function write_output_paper(data, sim_data, col_long, fp, param_vector, obj_eval, output_dir, paper_dir, data_mom, sim_mom)

% This file combines code from all other mfiles that generate latex code: descrip_stat, technology_fractions, policy_graphs. etc.
% In order to modify the tables accordingly, we generate a different
% txt-file for every table we want to include in the final version of the
% paper. We also write the necessary latex code to the screen here.

% Set the current directory to where the paper output is saved;

cd(paper_dir)

% Load vector of Bootstrapped parameters and SEs for tables in paper
filename = [output_dir,'temp_param_bootstrap.mat'];
if exist(filename,'file')==0
    disp('Vector of bootstrapped model parameter estimates not found!')
    disp('Need to rerun bootstrap first by setting fp.do_boot_se equal to 1 in master file.')
    disp('Assuming standard errors are all 0 for now so the code can compile.')
    vector_SE = zeros(150,1);
else % load existing vector of bootstrap parameters
    load([output_dir,'temp_param_bootstrap'],'param_est_boot_matrix');
    % calculate SEs for all transformed parameters
    vector_SE = nanstd(param_est_boot_matrix,1,2);
end

% Load SEs for the NLI process, which is estimated separately

filename2 = [output_dir,'temp_param_bootstrap_NLI.mat'];
if exist(filename2,'file')==0
    disp('Vector of bootstrapped non-labor income parameter estimates not found!')
    disp('Need to rerun bootstrap first by setting fp.do_boot_se and redo_BS_NLI equal to 1 in master file.')
    disp('Assuming standard errors are all 0 for now so the code can compile.')
    vector_SE_NLI = zeros(20,1);
else % load existing vector of bootstrap NLI parameters
    load([output_dir 'temp_param_bootstrap_NLI'],'param_est_boot_matrix_NLI');
    % calculate SEs for all transformed parameters
    vector_SE_NLI = nanstd(param_est_boot_matrix_NLI,1,2);
end

% Confidence intervals for NLI parameters:
% confidence_interval_SE_NLI = [prctile(param_est_boot_matrix_NLI,2.5,2) prctile(param_est_boot_matrix_NLI,5,2) ...
%     prctile(param_est_boot_matrix_NLI,95,2) prctile(param_est_boot_matrix_NLI,97.5,2)]; % 4 percentiles in increasing order

index97 = data(:,col_long.age) == data(:,col_long.child_age97);
index02 = data(:,col_long.age) == data(:,col_long.child_age97)+5;
index07 = data(:,col_long.age) == data(:,col_long.child_age97)+10;

%% Table: Descriptive Statistics

mom_age = data(index97,col_long.mom_age);
dad_age = data(index97,col_long.dad_age);
mom_educ = data(index97,col_long.mom_educ);
dad_educ = data(index97,col_long.dad_educ);
child_age97 = data(index97,col_long.age);

% child_age97 = data(index97,col_long.age);
LW97 = data(index97,col_long.lw_raw);
nrobs_LW97 = sum(not(isnan(LW97)));

LW02 = data(index02,col_long.lw_raw);
nrobs_LW02 = sum(not(isnan(LW02)));

LW07 = data(index07,col_long.lw_raw);
nrobs_LW07 = sum(not(isnan(LW07)));

h1 = data(:,col_long.mom_hours);
h2 = data(:,col_long.dad_hours);
w1 = data(:,col_long.mom_wage);
w2 = data(:,col_long.dad_wage);
NLI = data(:,col_long.nonlabor_income);

prec = '%.2f'; % precision for this table (except column "N")

diary off;
diary_dir = [cd,'\table_descrip_stat.txt'];
force_delete(diary_dir);
diary('table_descrip_stat.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Descriptive Statistics}']);
disp(['\label{table:descrip_stat}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lrrr}']);
disp(['\hline \hline']);
disp(['\multicolumn{4}{c}{\rule{0pt}{3ex} \textbf{PSID-CDS: 1997, 2002, 2007}} \\']);
disp(['\hline']);
disp(['\rule{0pt}{3ex} & \textbf{Mean} & \textbf{Std.} & \textbf{N} \\ ']);

disp(['Mother''s age in 1997 & ',num2str(mean(mom_age),prec),' & ',num2str(std(mom_age),prec), ' & ',num2str(sum(mom_age>0),'%.0f'),' \\', ...
    'Father''s age in 1997 & ',num2str(mean(dad_age),prec),' & ',num2str(std(dad_age),prec), ' & ',num2str(sum(dad_age>0),'%.0f'),' \\', ...
    'Mother''s education & ',num2str(mean(mom_educ),prec),' & ',num2str(std(mom_educ),prec), ' & ',num2str(sum(mom_educ>0),'%.0f'),' \\', ...
    'Father''s education & ',num2str(mean(dad_educ),prec),' & ',num2str(std(dad_educ),prec), ' & ',num2str(sum(dad_educ>0),'%.0f'),' \\', ...;
    'Child''s age in 1997 & ',num2str(mean(child_age97),prec),' & ',num2str(std(child_age97),prec), ' & ',num2str(sum(child_age97>0),'%.0f'),' \\', ...
    'Letter Word raw score in 1997 & ',num2str(mean(LW97),prec),' & ',num2str(std(LW97),prec), ' & ',num2str(nrobs_LW97,'%.0f'),' \\', ...
    'Letter Word raw score in 2002 & ',num2str(mean(LW02),prec),' & ',num2str(std(LW02),prec), ' & ',num2str(nrobs_LW02,'%.0f'),' \\', ...
    'Letter Word raw score in 2007 & ',num2str(nanmean(LW07),prec),' & ',num2str(nanstd(LW07),prec), ' & ',num2str(nrobs_LW07,'%.0f'),' \\']);

disp(['\hline']);
disp(['\multicolumn{4}{c}{\rule{0pt}{3ex} \textbf{PSID-Core: 1986 - 2010}} \\']);
disp(['\hline']);
disp(['\rule{0pt}{3ex} & \textbf{Mean} & \textbf{Std.} & \textbf{N} \\ ']);

disp(['Mother''s work hours per week & ',num2str(nanmean(h1),prec),' & ',num2str(nanstd(h1),prec), ' & ',num2str(sum(not(isnan(h1))),'%.0f'),' \\', ...
    'Father''s work hours per week & ',num2str(nanmean(h2),prec),' & ',num2str(nanstd(h2),prec), ' & ',num2str(sum(not(isnan(h2))),'%.0f'),' \\', ...
    'Mother''s hourly wage & ',num2str(nanmean(w1),prec),' & ',num2str(nanstd(w1),prec), ' & ',num2str(sum(not(isnan(w1))),'%.0f'),' \\', ...
    'Father''s hourly wage & ',num2str(nanmean(w2),prec),' & ',num2str(nanstd(w2),prec), ' & ',num2str(sum(not(isnan(w2))),'%.0f'),' \\', ...
    'Non-labor income per week & ',num2str(nanmean(NLI),prec),' & ',num2str(nanstd(NLI),prec), ' & ',num2str(sum(not(isnan(NLI))),'%.0f'),' \\']);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes}: Parental work hours, wages and non-labor income statistics are averaged over all years where the child is between 0 and 16 years old, ranging from 1986 to 2010.\\']);
disp(['\textbf{Source}: PSID-CDS combined sample from 1997, 2002 and 2007 interviews and PSID core data between 1986 and 2010.']);
disp(['\end{table}']);

diary off;

%% Table: Time Allocation by Child Age

% Consider 4 different age ranges
index1 = data(:,col_long.age) >= 3 & data(:,col_long.age) <= 5; % pre-school/kindergarten
index2 = data(:,col_long.age) >= 6 & data(:,col_long.age) <= 10; % elementary age
index3 = data(:,col_long.age) >= 11 & data(:,col_long.age) <= 13; % middle school age
index4 = data(:,col_long.age) >= 14 & data(:,col_long.age) <= 16; % high school age

h1 = data(:,col_long.mom_hours);
h2 = data(:,col_long.dad_hours);
tau1 = data(:,col_long.tau1);
tau2 = data(:,col_long.tau2);
tau_c = data(:,col_long.tau_c);
tau_school = data(:,col_long.tau_school);
% leis1 = fp.TT_p - tau1 - h1; % dad leisure in data
% leis2 = fp.TT_p - tau2 - h2; % mom leisure in data
% Note: parental leisure at ages 3-7 is NaN, because all those
% observations are from CDS-I in 1997, when we don't observe h1/h2 since
% there was no PSID interview in that year. Therefore, we define leis1 and
% leis2 as the residual.
% Child leisure:
leisc = fp.TT_p - tau_school-tau_c-tau1-tau2; % child leisure, using tau_school at the individual level

prec = '%.2f'; % precision for this table

diary off;
diary_dir = [cd,'\table_time_alloc.txt'];
force_delete(diary_dir);
diary('table_time_alloc.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Time Allocation by Child Age (Average Hours per Week)}']);
disp(['\label{table:time_alloc}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lrrrr}']);
disp(['\hline \hline']);
disp(['\textbf{Child Age} \rule{0pt}{3ex}& \multicolumn{1}{c}{\textbf{3-5}} & ',...
    '\multicolumn{1}{c}{\textbf{6-10}} & ',...
    '\multicolumn{1}{c}{\textbf{11-13}} & ',...
    '\multicolumn{1}{c}{\textbf{14-16}} \\']);
disp(['\hline']);

disp([
    'Mother''s Work Hours, Avg. \rule{0pt}{3ex}& ',...
    num2str(nanmean(h1(index1,:)),prec), ' & ',...
    num2str(nanmean(h1(index2,:)),prec), ' & ',...
    num2str(nanmean(h1(index3,:)),prec), ' & ', ...
    num2str(nanmean(h1(index4,:)),prec), ' \\ ',...    
    'Mother''s Work Hours, St. Dev. & ',...
    num2str(nanstd(h1(index1,:)),prec), ' & ',...
    num2str(nanstd(h1(index2,:)),prec), ' & ',...
    num2str(nanstd(h1(index3,:)),prec), ' & ', ...
    num2str(nanstd(h1(index4,:)),prec), ' \\ ',... 
    'Mother''s Investment Time, Avg. & ',...
    num2str(nanmean(tau1(index1,:)),prec), ' & ',...
    num2str(nanmean(tau1(index2,:)),prec), ' & ',...
    num2str(nanmean(tau1(index3,:)),prec), ' & ', ...
    num2str(nanmean(tau1(index4,:)),prec), ' \\ ',...
    'Mother''s Investment Time, St. Dev. & ',...
    num2str(nanstd(tau1(index1,:)),prec), ' & ',...
    num2str(nanstd(tau1(index2,:)),prec), ' & ',...
    num2str(nanstd(tau1(index3,:)),prec), ' & ', ...
    num2str(nanstd(tau1(index4,:)),prec), ' \\ ',...
    'Mother''s Leisure, Avg. (residual) & ',...
    num2str(fp.TT_p - nanmean(h1(index1,:))- nanmean(tau1(index1,:)),prec), ' & ',...
    num2str(fp.TT_p - nanmean(h1(index2,:))- nanmean(tau1(index2,:)),prec), ' & ',...
    num2str(fp.TT_p - nanmean(h1(index3,:))- nanmean(tau1(index3,:)),prec), ' & ',...
    num2str(fp.TT_p - nanmean(h1(index4,:))- nanmean(tau1(index4,:)),prec), ' \\ ',...
    'Father''s Work Hours, Avg. \rule{0pt}{3ex} & ',...
    num2str(nanmean(h2(index1,:)),prec), ' & ',...
    num2str(nanmean(h2(index2,:)),prec), ' & ',...
    num2str(nanmean(h2(index3,:)),prec), ' & ', ...
    num2str(nanmean(h2(index4,:)),prec), ' \\ ',...
    'Father''s Work Hours, St. Dev. & ',...
    num2str(nanstd(h2(index1,:)),prec), ' & ',...
    num2str(nanstd(h2(index2,:)),prec), ' & ',...
    num2str(nanstd(h2(index3,:)),prec), ' & ', ...
    num2str(nanstd(h2(index4,:)),prec), ' \\ ',...
    'Father''s Investment Time, Avg. & ',...
    num2str(nanmean(tau2(index1,:)),prec), ' & ',...
    num2str(nanmean(tau2(index2,:)),prec), ' & ',...
    num2str(nanmean(tau2(index3,:)),prec), ' & ', ...
    num2str(nanmean(tau2(index4,:)),prec), ' \\ ',...
    'Father''s Investment Time, St. Dev. & ',...
    num2str(nanstd(tau2(index1,:)),prec), ' & ',...
    num2str(nanstd(tau2(index2,:)),prec), ' & ',...
    num2str(nanstd(tau2(index3,:)),prec), ' & ', ...
    num2str(nanstd(tau2(index4,:)),prec), ' \\ ',...
    'Father''s Leisure, Avg. (residual) & ',...
    num2str(fp.TT_p - nanmean(h2(index1,:))- nanmean(tau2(index1,:)),prec), ' & ',...
    num2str(fp.TT_p - nanmean(h2(index2,:))- nanmean(tau2(index2,:)),prec), ' & ',...
    num2str(fp.TT_p - nanmean(h2(index3,:))- nanmean(tau2(index3,:)),prec), ' & ',...
    num2str(fp.TT_p - nanmean(h2(index4,:))- nanmean(tau2(index4,:)),prec), ' \\ ',...   
    'Child''s Self-Investment Time, Avg. \rule{0pt}{3ex} & ',...
    num2str(nanmean(tau_c(index1,:)),prec), ' & ',...
    num2str(nanmean(tau_c(index2,:)),prec), ' & ',...
    num2str(nanmean(tau_c(index3,:)),prec), ' & ', ...
    num2str(nanmean(tau_c(index4,:)),prec), ' \\ ',...
    'Child''s Self-Investment Time, St. Dev. & ',...
    num2str(nanstd(tau_c(index1,:)),prec), ' & ',...
    num2str(nanstd(tau_c(index2,:)),prec), ' & ',...
    num2str(nanstd(tau_c(index3,:)),prec), ' & ', ...
    num2str(nanstd(tau_c(index4,:)),prec), ' \\ ',...
    'Child''s School Time, Avg. & ',...
    num2str(nanmean(tau_school(index1,:)),prec), ' & ',...
    num2str(nanmean(tau_school(index2,:)),prec), ' & ',...
    num2str(nanmean(tau_school(index3,:)),prec), ' & ', ...
    num2str(nanmean(tau_school(index4,:)),prec), ' \\ ',...
    'Child''s Leisure, Avg. (residual) & ',...
    num2str(nanmean(leisc(index1,:)),prec), ' & ',...
    num2str(nanmean(leisc(index2,:)),prec), ' & ',...
    num2str(nanmean(leisc(index3,:)),prec), ' & ', ...
    num2str(nanmean(leisc(index4,:)),prec), ' \\ ']);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes}: All time use statistics are calculated over all years where the child is between 3 and 16 years old. ',...
    'Parental leisure is defined as an average residual, after subtracting average work hours and investment time from the total weekly time budget of $112$ hours. ',...
    'Child leisure is defined as a residual, after subtracting parental investment time, self-investment time, and school time from the total weekly time budget of $112$ hours. \\']);
disp(['\textbf{Source}: PSID-CDS combined sample from 1997, 2002 and 2007 interviews and PSID core data between 1986 and 2010.']);
disp(['\end{table}']);

diary off;

%% Table: Parental labor supply by child age

% Define child age categories (proxy for preschool/elementary/middle/high)
beg_age = [3 6 11 14];
end_age = [5 10 13 16];

tmp_table = NaN(length(beg_age),20);

for t = 1:length(beg_age)
    % select the rows we need for this age group
    row_ind = (data(:,col_long.age) >= beg_age(t) & data(:,col_long.age) <= end_age(t));
    
    % LFP mothers
    temp1 = nansum(data( row_ind,col_long.mom_hours) > 0);
    temp2 = nansum((data( row_ind,col_long.mom_hours)) == 0);
    tmp_table(t,1) = temp1./(temp1 + temp2);
    
    % LFP fathers
    temp1 = nansum(data( row_ind,col_long.dad_hours) > 0);
    temp2 = nansum((data( row_ind,col_long.dad_hours)) == 0);
    tmp_table(t,2) = temp1./(temp1 + temp2);
    
    % Joint LFP
    temp1 = nansum( data( row_ind,col_long.dad_hours) > 0 & data(row_ind,col_long.mom_hours) > 0 );
    temp2 = nansum( data( row_ind,col_long.dad_hours) == 0 | data(row_ind,col_long.mom_hours) == 0 );
    tmp_table(t,3) = temp1./(temp1 + temp2);
    
    temp = (row_ind == 1 & data(:,col_long.mom_hours) > 0); % conditional on supplying positive hours
    tmp_table(t,4) = nanmean(data(temp,col_long.mom_hours));
    tmp_table(t,5) = nanstd(data(temp,col_long.mom_hours));
    
    temp = (row_ind == 1 & data(:,col_long.dad_hours) > 0); % conditional on supplying positive hours
    tmp_table(t,6) = nanmean(data(temp,col_long.dad_hours));
    tmp_table(t,7) = nanstd(data(temp,col_long.dad_hours));
end

prec1 = '%.3f'; % precision for  LFP
prec = '%.2f'; % precision for bottom half

diary off;
diary_dir = [cd,'\table_labor.txt'];
force_delete(diary_dir);
diary('table_labor.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Parental Labor Supply by Child Age}']);
disp(['\label{table:labor}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcccc}']);
disp(['\hline \hline']);
disp(['\textbf{Child Age} \rule{0pt}{3ex}& \multicolumn{1}{c}{\textbf{3-5}} & ',...
    '\multicolumn{1}{c}{\textbf{6-10}} & ',...
    '\multicolumn{1}{c}{\textbf{11-13}} & ',...
    '\multicolumn{1}{c}{\textbf{14-16}} \\']);
disp(['\hline']);
disp([' & \multicolumn{4}{c}{\rule{0pt}{3ex}\textbf{Fraction Working $> 0$ Hours}} \\']);

disp(['Mothers working \rule{0pt}{3ex}& ',num2str(tmp_table(1,1),prec1), ' & ',...
    num2str(tmp_table(2,1),prec1), ' & ',...
    num2str(tmp_table(3,1),prec1), ' & ', ...
    num2str(tmp_table(4,1),prec1), ' \\ ',...
    'Fathers working & ',num2str(tmp_table(1,2),prec1), ' & ',...
    num2str(tmp_table(2,2),prec1), ' & ',...
    num2str(tmp_table(3,2),prec1), ' & ', ...
    num2str(tmp_table(4,2),prec1), ' \\ ',...
    'Both working & ',num2str(tmp_table(1,3),prec1), ' & ',...
    num2str(tmp_table(2,3),prec1), ' & ',...
    num2str(tmp_table(3,3),prec1), ' & ', ...
    num2str(tmp_table(4,3),prec1), ' \\ ']);

disp([' & \multicolumn{4}{c}{\rule{0pt}{3ex}\textbf{Working Hours (cond. on Hours $> 0$)}} \\']);

disp(['Mothers, Avg. \rule{0pt}{3ex}& ',num2str(tmp_table(1,4),prec), ' & ',...
    num2str(tmp_table(2,4),prec), ' & ',...
    num2str(tmp_table(3,4),prec), ' & ', ...
    num2str(tmp_table(4,4),prec), ' \\ ',...
    'Mothers, St. Dev. & ',num2str(tmp_table(1,5),prec), ' & ',...
    num2str(tmp_table(2,5),prec), ' & ',...
    num2str(tmp_table(3,5),prec), ' & ', ...
    num2str(tmp_table(4,5),prec), ' \\ ',...
    'Fathers, Avg. & ',num2str(tmp_table(1,6),prec), ' & ',...
    num2str(tmp_table(2,6),prec), ' & ',...
    num2str(tmp_table(3,6),prec), ' & ', ...
    num2str(tmp_table(4,6),prec), ' \\ ',...
    'Fathers, St. Dev. & ',num2str(tmp_table(1,7),prec), ' & ',...
    num2str(tmp_table(2,7),prec), ' & ',...
    num2str(tmp_table(3,7),prec), ' & ', ...
    num2str(tmp_table(4,7),prec), ' \\ ']);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes}: Upper half of the table shows labor force participation rates. Bottom half shows averages and standard deviations of labor hours conditional on working positive hours.\\']);
disp(['\textbf{Source}: PSID-CDS combined sample from 1997, 2002 and 2007 interviews and PSID core data between 1986 and 2010.']);
disp(['\end{table}']);

diary off;

%% Summary statistics on (Un)conditional Allowance data: Joint probability of "giving allowance contingent on school work" for multiple different age categories

% define child age categories to reflect preschool/elementary/middle/high:
beg_age = [8 11 14 1];
end_age = [10 13 16 16];

beg_educ = [1 1 13 17];
end_educ = [17 12 16 17];
% Recall: allowance_ABCD2 means parental reports (more complete, reported for ages 6+) and
% allowance_ABCD means child reports (only reported for ages 12+)

tmp_table = NaN(length(beg_educ)+1,length(beg_age),3);

tmp_CCT = NaN(size(data,1),1);
% only look at 2002 + 2007 data when we observe conditional allowances
tmp0a = data(:,col_long.allowance_A2) == 0 & data(:,col_long.year) > 97; % question D might be missing here
tmp0b = data(:,col_long.allowance_A2) == 1 & data(:,col_long.allowance_D2) == 0 & data(:,col_long.year) > 97;
tmp1 = data(:,col_long.allowance_A2) == 1 & data(:,col_long.allowance_D2) == 1 & data(:,col_long.year) > 97;
% this excludes observations where A = 1 (i.e. they use allowance) but D is missing!
tmp_CCT(tmp0a==1)=0;
tmp_CCT(tmp0b==1)=0;
tmp_CCT(tmp1==1)=1;

for t = 1:length(beg_age)
    for jj = 1:length(beg_educ)
        % select the rows we need for this age group
        row_ind = (data(:,col_long.age) >= beg_age(t) & data(:,col_long.age) <= end_age(t) & data(:,col_long.dad_educ) >= beg_educ(jj) & data(:,col_long.dad_educ) <= end_educ(jj) & data(:,col_long.year) > 97);
        % Only include 2002 and 2007 data because we don't observe question D) in 1997!
        
        % A) Does your child get an allowance? Yes or no
        temp1 = nansum(data( row_ind,col_long.allowance_A2) == 1);
        temp2 = nansum(data( row_ind,col_long.allowance_A2) == 0);
        tmp_table(jj,t,1) = temp1/(temp1+temp2);
        
        % D) Is the allowance contingent on child doing his/her school work? Yes or no (only answered for 2002 and 2007)
        temp1 = nansum(data( row_ind,col_long.allowance_D2) == 1);
        temp2 = nansum(data( row_ind,col_long.allowance_D2) == 0);
        tmp_table(jj,t,2) = temp1/(temp1+temp2);
        
        % Joint probability: get an allowance AND it's contingent on school work (this drops 1997 data since question D wasn't asked then)
        tmp_table(jj,t,3) = nanmean(tmp_CCT(row_ind));
        if jj == 1 % fill out sample sizes whenever we look at the full educ range
            % A) % sample size
            tmp_table(length(beg_educ)+1,t,1) = sum(~isnan(data( row_ind,col_long.allowance_A2)));
            % D) sample size
            tmp_table(length(beg_educ)+1,t,2) = sum(~isnan(data( row_ind,col_long.allowance_D2)));
            % Joint sample size
            tmp_table(length(beg_educ)+1,t,3) = sum(~isnan(tmp_CCT(row_ind)));
        end
    end
end

% disp('Allowance use by child age and parental educ/income:')
% tmp_table
% tmp_table(1:4,:,1).*tmp_table(1:4,:,2) % joint probabilities

% Create table for paper

prec1 = '%.3f'; % precision
prec0 = '%.0f'; % precision for bottom half

diary off;
diary_dir = [cd,'\table_allowances2.txt'];
force_delete(diary_dir);
diary('table_allowances2.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\captionsetup{width=\linewidth}']);
disp(['\caption{(Un)conditional Allowances, by Child Age and Household Characteristics}']);
disp(['\label{table:allowances}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcccc}']);
disp(['\hline \hline']);
disp([' & \multicolumn{4}{c}{\rule{0pt}{3ex}\textbf{Fraction }} \\']);
disp(['\textbf{Child Age} & ',...
    '\multicolumn{1}{c}{\textbf{8-10}} & ',...
    '\multicolumn{1}{c}{\textbf{11-13}} & ',...
    '\multicolumn{1}{c}{\textbf{14-16}} & ',...
    '\multicolumn{1}{c}{\textbf{All ages}} \\']);
disp(['\hline']);

disp(['\multicolumn{5}{l}{\textbf{(a) Does your child receive an allowance?}}  \rule{0pt}{3ex}\\']);

disp(['All households \rule{0pt}{3ex}& ',...
    num2str(tmp_table(1,1,1),prec1), ' & ', num2str(tmp_table(1,2,1),prec1), ' & ', num2str(tmp_table(1,3,1),prec1), ' & ', num2str(tmp_table(1,4,1),prec1), ' \\ ',...
    'High school or less \rule{0pt}{3ex}& ',...
    num2str(tmp_table(2,1,1),prec1), ' & ', num2str(tmp_table(2,2,1),prec1), ' & ', num2str(tmp_table(2,3,1),prec1), ' & ', num2str(tmp_table(2,4,1),prec1), ' \\ ',...
    '(Some) college & ',...
    num2str(tmp_table(3,1,1),prec1), ' & ', num2str(tmp_table(3,2,1),prec1), ' & ', num2str(tmp_table(3,3,1),prec1), ' & ', num2str(tmp_table(3,4,1),prec1), ' \\ '...
    'Graduate & ',...
    num2str(tmp_table(4,1,1),prec1), ' & ', num2str(tmp_table(4,2,1),prec1), ' & ', num2str(tmp_table(4,3,1),prec1), ' & ', num2str(tmp_table(4,4,1),prec1), ' \\ '...
    'Sample size \rule{0pt}{3ex}& ', num2str(tmp_table(end,1,1),prec0), ' & ', num2str(tmp_table(end,2,1),prec0), ' & ', num2str(tmp_table(end,3,1),prec0), ' & ', num2str(tmp_table(end,4,1),prec0), ' \\ ']);

disp(['\multicolumn{5}{l}{\textbf{(b) Is the allowance contingent on child doing his/her school work?}}  \rule{0pt}{4ex}\\']);

disp(['All households \rule{0pt}{3ex}& ',...
    num2str(tmp_table(1,1,2),prec1), ' & ', num2str(tmp_table(1,2,2),prec1), ' & ', num2str(tmp_table(1,3,2),prec1), ' & ', num2str(tmp_table(1,4,2),prec1), ' \\ ',...
    'High school or less \rule{0pt}{3ex}& ',...
    num2str(tmp_table(2,1,2),prec1), ' & ', num2str(tmp_table(2,2,2),prec1), ' & ', num2str(tmp_table(2,3,2),prec1), ' & ', num2str(tmp_table(2,4,2),prec1), ' \\ ',...
    '(Some) college & ',...
    num2str(tmp_table(3,1,2),prec1), ' & ', num2str(tmp_table(3,2,2),prec1), ' & ', num2str(tmp_table(3,3,2),prec1), ' & ', num2str(tmp_table(3,4,2),prec1), ' \\ '...
    'Graduate & ',...
    num2str(tmp_table(4,1,2),prec1), ' & ', num2str(tmp_table(4,2,2),prec1), ' & ', num2str(tmp_table(4,3,2),prec1), ' & ', num2str(tmp_table(4,4,2),prec1), ' \\ '...
    'Sample size \rule{0pt}{3ex}& ', num2str(tmp_table(end,1,2),prec0), ' & ', num2str(tmp_table(end,2,2),prec0), ' & ', num2str(tmp_table(end,3,2),prec0), ' & ', num2str(tmp_table(end,4,2),prec0), ' \\ ']);

disp(['\multicolumn{5}{l}{\textbf{(c) Joint probability of giving an allowance conditional on school work}}  \rule{0pt}{4ex}\\']);

disp(['All households \rule{0pt}{3ex}& ',...
    num2str(tmp_table(1,1,3),prec1), ' & ', num2str(tmp_table(1,2,3),prec1), ' & ', num2str(tmp_table(1,3,3),prec1), ' & ', num2str(tmp_table(1,4,3),prec1), ' \\ ',...
    'High school or less \rule{0pt}{3ex}& ',...
    num2str(tmp_table(2,1,3),prec1), ' & ', num2str(tmp_table(2,2,3),prec1), ' & ', num2str(tmp_table(2,3,3),prec1), ' & ', num2str(tmp_table(2,4,3),prec1), ' \\ ',...
    '(Some) college & ',...
    num2str(tmp_table(3,1,3),prec1), ' & ', num2str(tmp_table(3,2,3),prec1), ' & ', num2str(tmp_table(3,3,3),prec1), ' & ', num2str(tmp_table(3,4,3),prec1), ' \\ '...
    'Graduate & ',...
    num2str(tmp_table(4,1,3),prec1), ' & ', num2str(tmp_table(4,2,3),prec1), ' & ', num2str(tmp_table(4,3,3),prec1), ' & ', num2str(tmp_table(4,4,3),prec1), ' \\ '...
    'Sample size \rule{0pt}{3ex}& ', num2str(tmp_table(end,1,3),prec0), ' & ', num2str(tmp_table(end,2,3),prec0), ' & ', num2str(tmp_table(end,3,3),prec0), ' & ', num2str(tmp_table(end,4,3),prec0), ' \\ ']);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes}: The table shows the fraction of households who answer ``Yes" to the questions: ',...
    '(a) ``Does your child receive an allowance?", and (b) ``Is the allowance contingent on your child doing his/her school work?". ',...
    'Panel (c) shows the joint probabilities of using a conditional allowance, defined as the product of ',...
    'the corresponding fractions in panel (a) and (b). The various rows break down the responses by the father''s ',...
    'education level (high school or less, some college or college degree, or graduate level). ',...
    'Note: Question (b) was not asked in 1997, and therefore all numbers only include data from 2002 and 2007, when all children were at least 8 years old.\\']);
disp('\textbf{Source}: PSID-CDS combined sample from 2002 and 2007.');
disp(['\end{table}']);

diary off;

%% Table: Preference Parameters

prec1 = '%.3f'; % precision for estimates
prec2 = '%.5f'; % precision for standard errors

mean_alpha1 = param_vector(1);
mean_alpha2 = param_vector(2);
mean_alpha3 = param_vector(3);
mean_alpha4 = param_vector(4);
mean_lambda1 = param_vector(5);
mean_lambda2 = param_vector(6);
mean_lambda3 = param_vector(7);

mean_alpha1_SE = vector_SE(1);
mean_alpha2_SE = vector_SE(2);
mean_alpha3_SE = vector_SE(3);
mean_alpha4_SE = vector_SE(4);
mean_lambda1_SE = vector_SE(5);
mean_lambda2_SE = vector_SE(6);
mean_lambda3_SE = vector_SE(7);
    
phi = param_vector(39);
cost_mean = param_vector(89);
cost_mean_SE = vector_SE(89);

diary off;
diary_dir = [cd,'\table_param_prefs.txt'];
force_delete(diary_dir);
diary('table_param_prefs.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Preference Parameter Estimates}']);
disp(['\label{table:param_prefs}']);

disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcc}']);
disp(['\hline \hline']);

disp(['\multicolumn{1}{l}{\textbf{(a) Parental preferences}}  \rule{0pt}{3ex}& \multicolumn{1}{c}{\textbf{Estimate}} & \multicolumn{1}{c}{\textbf{SE}}\\']);
disp(['Mother''s Leisure ($\alpha_1$) & ', num2str(mean_alpha1,prec1), ' & (', num2str(mean_alpha1_SE,prec2), ') \\ ',...
    'Father''s Leisure ($\alpha_2$) & ', num2str(mean_alpha2,prec1), ' & (', num2str(mean_alpha2_SE,prec2), ') \\ ',...
    'Consumption ($\alpha_3$) & ', num2str(mean_alpha3,prec1), ' & (', num2str(mean_alpha3_SE,prec2), ') \\ ',...
    'Child Quality ($\alpha_4$) & ', num2str(mean_alpha4,prec1), ' & (', num2str(mean_alpha4_SE,prec2), ') \\ ']);

disp(['\multicolumn{1}{l}{\textbf{(b) Child preferences}}  \rule{0pt}{4ex}& \multicolumn{1}{c}{\textbf{Estimate}} & \multicolumn{1}{c}{\textbf{SE}}\\']);
disp(['Leisure ($\lambda_1$) & ', num2str(mean_lambda1,prec1), ' & (', num2str(mean_lambda1_SE,prec2), ') \\ ',...
    'Consumption ($\lambda_2$) & ', num2str(mean_lambda2,prec1), ' & (', num2str(mean_lambda2_SE,prec2), ') \\ ',...
    'Child Quality ($\lambda_3$) & ', num2str(mean_lambda3,prec1), ' & (', num2str(mean_lambda3_SE,prec2), ') \\ ']);

disp('\multicolumn{1}{l}{\textbf{(c) Other preferences}} \rule{0pt}{4ex}& \multicolumn{1}{c}{\textbf{Estimate}} & \multicolumn{1}{c}{\textbf{SE}}  \\');
disp(['Parental altruism ($\varphi$) & ',num2str(phi,prec1), ' & ', '-', ' \\ ',...
    'Mean/Std. of CCT Utility Cost ($\kappa$) & ',num2str(cost_mean,prec1),' & (', num2str(cost_mean_SE,prec2), ') \\ ']);

disp('\hline \hline');
disp('\end{tabular*}');
disp('\end{center}');
disp(['\textbf{Notes}: SEs are standard errors computed using a cluster bootstrap ',...
    'sampling each household with replacement. Parameters without SE are assumed (not estimated) values.']);
disp('\end{table}');

diary off;

%% Table: Cognitive/Noncognitive Technology Parameters

prec1 = '%.3f'; % precision for estimates
prec2 = '%.5f'; % precision for standard errors

delta1_intercept = param_vector(40);
delta1_slope = param_vector(41);
delta1_educ = param_vector(84); % tau1, mom_educ
delta2_intercept = param_vector(42);
delta2_slope = param_vector(43);
delta2_educ = param_vector(85); % tau2, dad_educ
delta3_intercept = param_vector(44);
delta3_slope = param_vector(45);
delta4_intercept = param_vector(46);
delta4_slope = param_vector(47);
delta5_intercept = param_vector(48);
delta5_slope = param_vector(49);

%SEs
delta1_intercept_SE = vector_SE(40);
delta1_slope_SE = vector_SE(41);
delta1_educ_SE = vector_SE(84); % tau1, mom_educ
delta2_intercept_SE = vector_SE(42);
delta2_slope_SE = vector_SE(43);
delta2_educ_SE = vector_SE(85); % tau2, dad_educ
delta3_intercept_SE = vector_SE(44);
delta3_slope_SE = vector_SE(45);
delta4_intercept_SE = vector_SE(46);
delta4_slope_SE = vector_SE(47);
delta5_intercept_SE = vector_SE(48);
delta5_slope_SE = vector_SE(49);


% Load TFP parameters
% TFP_vec = [1;1;param_vector(90:103)]; % age 1-2: TFP = 1, age 3-16: 14 years
% TFP_vec_SE = [0;0;vector_SE(90:103)];

TFP_vec2 = param_vector(104:107); % L0,L,k,t0
TFP_vec2_SE = vector_SE(104:107);

% Load parameters characterizing the Markov transition matrix
b0_up = param_vector(108);
b1_up = param_vector(109);
b2_up = param_vector(110);
b3_up = param_vector(111);

b0_down = param_vector(112);
b1_down = param_vector(113);
b2_down = param_vector(114); % assumed symmetry
b3_down = param_vector(115); % assumed symmetry

b0_up_SE = vector_SE(108);
b1_up_SE = vector_SE(109);
b2_up_SE = vector_SE(110);
b3_up_SE = vector_SE(111);

b0_down_SE = vector_SE(112);
b1_down_SE = vector_SE(113);
%     b2_down = param_vector(114); % assumed symmetry - not estimated
%     b3_down = param_vector(115); % assumed symmetry - not estimated

% Load parameters characterizing the initial conditions for \beta_{c,t0}
nu10 = param_vector(116);
nu11 = param_vector(117);
nu12 = param_vector(118);
nu20 = param_vector(119);
nu21 = param_vector(120);
nu22 = param_vector(121);

nu10_SE = vector_SE(116);
nu11_SE = vector_SE(117);
nu12_SE = vector_SE(118);
nu20_SE = vector_SE(119);
nu21_SE = vector_SE(120);
nu22_SE = vector_SE(121);

diary off;
diary_dir = [cd,'\table_param_tech.txt'];
force_delete(diary_dir);
diary('table_param_tech.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Technology Parameter Estimates}']);
disp(['\label{table:param_tech}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} llcc}']);
disp(['\hline \hline']);

disp(['\multicolumn{4}{l}{\textbf{(a) Cognitive Human Capital}}  \rule{0pt}{3ex} \\ ']);
disp([' & & \multicolumn{1}{c}{\textbf{Estimate}} \rule{0pt}{3ex}& \multicolumn{1}{c}{\textbf{SE}}\\']);

disp(['Mother''s Active Time ($\delta_1$) & Intercept $d_{1,0}$&', num2str(delta1_intercept,prec1), ' & (', num2str(delta1_intercept_SE,prec2), ') \\',...
    ' & Slope - Age $d_{1,1}$ &', num2str(delta1_slope,prec1), ' & (', num2str(delta1_slope_SE,prec2), ') \\',...
    ' & Slope - Educ. $d_{1,2}$ &', num2str(delta1_educ,prec1), ' & (', num2str(delta1_educ_SE,prec2), ') \\',...
    'Father''s Active Time ($\delta_2$) & Intercept $d_{2,0}$ &', num2str(delta2_intercept,prec1), ' & (', num2str(delta2_intercept_SE,prec2), ') \\',...
    ' & Slope - Age $d_{2,1}$ &', num2str(delta2_slope,prec1), ' & (', num2str(delta2_slope_SE,prec2), ') \\',...
    ' & Slope - Educ. $d_{2,2}$ &', num2str(delta2_educ,prec1), ' & (', num2str(delta2_educ_SE,prec2), ') \\',...
    'Child''s Self-Investment Time ($\delta_3$) & Intercept $d_{3,0}$ &', num2str(delta3_intercept,prec1), ' & (', num2str(delta3_intercept_SE,prec2), ') \\',...
    ' & Slope - Age $d_{3,1}$ &', num2str(delta3_slope,prec1), ' & (', num2str(delta3_slope_SE,prec2), ') \\']);

disp(['Child Expenditures ($\delta_4$) & Intercept $d_{4,0}$ &', num2str(delta4_intercept,prec1), ' & (', num2str(delta4_intercept_SE,prec2), ') \\',...
    ' & Slope - Age $d_{4,1}$ &', num2str(delta4_slope,prec1), ' & (', num2str(delta4_slope_SE,prec2), ') \\',...
    'Last Period''s Child Quality ($\delta_5$) & Intercept $d_{5,0}$ &', num2str(delta5_intercept,prec1), ' & (', num2str(delta5_intercept_SE,prec2), ') \\',...
    ' & Slope - Age $d_{5,1}$ &', num2str(delta5_slope,prec1), ' & (', num2str(delta5_slope_SE,prec2), ') \\ ',...
    'Total Factor Productivity ($R_t$) \rule{0pt}{3ex} &  $d_{6,0}$', ' & ', num2str(TFP_vec2(1),prec2), ' & (', num2str(TFP_vec2_SE(1),prec2), ') \\',...
    '&  $d_{6,1}$', ' & ', num2str(TFP_vec2(2),prec2), ' & (', num2str(TFP_vec2_SE(2),prec2), ') \\',...
    '&  $d_{6,2}$', ' & ', num2str(TFP_vec2(3),prec2), ' & (', num2str(TFP_vec2_SE(3),prec2), ') \\',...
    '&  $d_{6,3}$', ' & ', num2str(TFP_vec2(4),prec2), ' & (', num2str(TFP_vec2_SE(4),prec2), ') \\ \\']);

disp(['\multicolumn{4}{l}{\textbf{(b) Child Discount Factors}}  \rule{0pt}{3ex} \\']);

disp(['Upward transitions ($D^{up}$) & Intercept $b^{up}_0$ &', num2str(b0_up,prec1), ' & (', num2str(b0_up_SE,prec2), ') \\',...
    ' & Slope - Age $b^{up}_1$ &', num2str(b1_up,prec1), ' & (', num2str(b1_up_SE,prec2), ') \\',...
    ' & Slope - CCT $b^{up}_2$ &', num2str(b2_up,prec1), ' & (', num2str(b2_up_SE,prec2), ') \\',...
    ' & Slope - Interaction $b^{up}_3$ &', num2str(b3_up,prec1), ' & (', num2str(b3_up_SE,prec2), ') \\',...
    'Downward transitions ($D^{down}$) & Intercept $b^{down}_0$ &', num2str(b0_down,prec1), ' & (', num2str(b0_down_SE,prec2), ') \\',...
    ' & Slope - Age $b^{down}_1$ &', num2str(b1_down,prec1), ' & (', num2str(b1_down_SE,prec2), ') \\',...
    ' & Slope - CCT $b^{down}_2$ &', num2str(b2_down,prec1), ' & - \\',...
    ' & Slope - Interaction $b^{down}_3$ &', num2str(b3_down,prec1), ' & - \\',...
    'Initial conditions - $Pr(\beta_{c,t^h_0} = \beta_c^1|t^h_0,s^h)$ & Intercept $\nu^1_0$ &', num2str(nu10,prec1), ' & (', num2str(nu10_SE,prec2), ') \\',...
    ' & Slope - Age $\nu^1_1$ &', num2str(nu11,prec1), ' & (', num2str(nu11_SE,prec2), ') \\',...
    ' & Slope - Parent Educ. $\nu^1_2$ &', num2str(nu12,prec1), ' & (', num2str(nu12_SE,prec2), ') \\',...
    'Initial conditions - $Pr(\beta_{c,t^h_0} = \beta_c^2|t^h_0,s^h)$ & Intercept $\nu^2_0$ &', num2str(nu20,prec1), ' & (', num2str(nu20_SE,prec2), ') \\',...
    ' & Slope - Age $\nu^2_1$ &', num2str(nu21,prec1), ' & (', num2str(nu21_SE,prec2), ') \\',...
    ' & Slope - Parent Educ. $\nu^2_2$ &', num2str(nu22,prec1), ' & (', num2str(nu22_SE,prec2), ') \\ \\']);

disp('\hline \hline');
disp('\end{tabular*}');
disp('\end{center}');

disp(['Notes: SEs are standard errors computed using a cluster bootstrap sampling each household with replacement.']);
disp('\end{table}');

diary off;

%% Table: Wage and NLI Parameters

prec1 = '%.3f'; % precision for estimates
prec2 = '%.5f'; % precision for standard errors

% Wage process parameters
beta0_mom = param_vector(54);
beta1_mom = param_vector(55);
beta2_mom = param_vector(56);
sigma_mom = param_vector(57);
beta0_dad = param_vector(58);
beta1_dad = param_vector(59);
beta2_dad = param_vector(60);
sigma_dad = param_vector(61);
wage_corr = param_vector(62);
% SEs
beta0_mom_SE = vector_SE(54);
beta1_mom_SE = vector_SE(55);
beta2_mom_SE = vector_SE(56);
sigma_mom_SE = vector_SE(57);
beta0_dad_SE = vector_SE(58);
beta1_dad_SE = vector_SE(59);
beta2_dad_SE = vector_SE(60);
sigma_dad_SE = vector_SE(61);
wage_corr_SE = vector_SE(62);

% Non-labor income parameters (estimated outside the model)
tmp = fp.NLI_par; % load optimal parameters

% step 1: logit to estimate Pr(I>0)
mu1 = tmp(1); % intercept
mu2 = tmp(2); % mom_age
mu3 = tmp(3); % mom_age_squared
mu4 = tmp(4); % mom_educ
mu5 = tmp(5); % dad_age
mu6 = tmp(6); % dad_age_squared
mu7 = tmp(7); % dad_educ
% step 2: NLI conditional on being positive
mu8 = tmp(8); % intercept
mu9 = tmp(9); % mom_age
mu10 = tmp(10); % mom_age_squared
mu11 = tmp(11); % mom_educ
mu12 = tmp(12); % dad_age
mu13 = tmp(13); % dad_age_squared
mu14 = tmp(14); % dad_educ
sigma_NLI = tmp(15); % st.dev. of epsilons

% step 1: logit to estimate Pr(I>0)
mu1_SE = vector_SE_NLI(1); % intercept
mu2_SE = vector_SE_NLI(2); % mom_age
mu3_SE = vector_SE_NLI(3); % mom_age_squared
mu4_SE = vector_SE_NLI(4); % mom_educ
mu5_SE = vector_SE_NLI(5); % dad_age
mu6_SE = vector_SE_NLI(6); % dad_age_squared
mu7_SE = vector_SE_NLI(7); % dad_educ

% step 2: NLI conditional on being positive
mu8_SE = vector_SE_NLI(8); % intercept
mu9_SE = vector_SE_NLI(9); % mom_age
mu10_SE = vector_SE_NLI(10); % mom_age_squared
mu11_SE = vector_SE_NLI(11); % mom_educ
mu12_SE = vector_SE_NLI(12); % dad_age
mu13_SE = vector_SE_NLI(13); % dad_age_squared
mu14_SE = vector_SE_NLI(14); % dad_educ
sigma_NLI_SE = vector_SE_NLI(15); % st.dev. of epsilons

diary off;
diary_dir = [cd,'\table_param_wages.txt'];
force_delete(diary_dir);
diary('table_param_wages.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Wage and Income Parameter Estimates}']);
disp(['\label{table:param_wages}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcccc}']);
disp(['\hline \hline']);

disp(['\multicolumn{1}{l}{\textbf{(a) Mother''s Log Wage Offer}} \rule{0pt}{3ex}& \multicolumn{1}{c}{\textbf{Estimate}} & \multicolumn{1}{c}{\textbf{SE}}\\']);

disp(['Intercept ($\eta_{0,1}$) & ', num2str(beta0_mom,prec1), ' & (', num2str(beta0_mom_SE,prec2), ') \\',...
    'Mother''s Age ($\eta_{1,1}$) & ', num2str(beta1_mom,prec1), ' & (', num2str(beta1_mom_SE,prec2), ') \\',...
    'Mother''s Education ($\eta_{2,1}$) & ', num2str(beta2_mom,prec1), ' & (', num2str(beta2_mom_SE,prec2), ') \\',...
    'Standard Deviation of Shock ($\sigma_{w_1}$) & ', num2str(sigma_mom,prec1), ' & (', num2str(sigma_mom_SE,prec2), ') \\', ...
    'Correlation with Father''s Wage Shock ($\rho_{12}$) & ', num2str(wage_corr,prec1), ' & (', num2str(wage_corr_SE,prec2), ') \\ ']);

disp(['\multicolumn{1}{l}{\textbf{(b) Father''s Log Wage Offer}} \rule{0pt}{4ex}& \multicolumn{1}{c}{\textbf{Estimate}} & \multicolumn{1}{c}{\textbf{SE}}\\']);

disp(['Intercept ($\eta_{0,2}$) & ', num2str(beta0_dad,prec1), ' & (', num2str(beta0_dad_SE,prec2), ') \\',...
    'Father''s Age ($\eta_{1,2}$) & ', num2str(beta1_dad,prec1), ' & (', num2str(beta1_dad_SE,prec2), ') \\',...
    'Father''s Education ($\eta_{2,2}$) & ', num2str(beta2_dad,prec1), ' & (', num2str(beta2_dad_SE,prec2), ') \\',...
    'Standard Deviation of Shock ($\sigma_{w_2}$) & ', num2str(sigma_dad,prec1), ' & (', num2str(sigma_dad_SE,prec2), ') \\ ']);

disp(['\multicolumn{1}{l}{\textbf{(c) Latent Non-Labor Income}} \rule{0pt}{4ex}& \multicolumn{1}{c}{\textbf{Estimate}} & \multicolumn{1}{c}{\textbf{SE}}\\']);

disp(['Logit - Intercept ($\mu_1$) & ', num2str(mu1,prec1), ' & (', num2str(mu1_SE,prec2), ') \\',...
    'Logit - Mother''s Age ($\mu_2$) & ', num2str(mu2,prec1), ' & (', num2str(mu2_SE,prec2), ') \\',...
    'Logit - Mother''s Age Squared ($\mu_3$) & ', num2str(mu3,prec1), ' & (', num2str(mu3_SE,prec2), ') \\',...
    'Logit - Mother''s Education ($\mu_4$) & ', num2str(mu4,prec1), ' & (', num2str(mu4_SE,prec2), ') \\',...
    'Logit - Father''s Age ($\mu_5$) & ', num2str(mu5,prec1), ' & (', num2str(mu5_SE,prec2), ') \\',...
    'Logit - Father''s Age Squared ($\mu_6$) & ', num2str(mu6,prec1), ' & (', num2str(mu6_SE,prec2), ') \\',...
    'Logit - Father''s Education ($\mu_7$) & ', num2str(mu7,prec1), ' & (', num2str(mu7_SE,prec2), ') \\',...
    'Conditional - Intercept ($\mu_8$) & ', num2str(mu8,prec1), ' & (', num2str(mu8_SE,prec2), ') \\',...
    'Conditional - Mother''s Age ($\mu_9$) & ', num2str(mu9,prec1), ' & (', num2str(mu9_SE,prec2), ') \\',...
    'Conditional - Mother''s Age Squared ($\mu_{10}$) & ', num2str(mu10,prec1), ' & (', num2str(mu10_SE,prec2), ') \\',...
    'Conditional - Mother''s Education ($\mu_{11}$) & ', num2str(mu11,prec1), ' & (', num2str(mu11_SE,prec2), ') \\',...
    'Conditional - Father''s Age ($\mu_{12}$) & ', num2str(mu12,prec1), ' & (', num2str(mu12_SE,prec2), ') \\',...
    'Conditional - Father''s Age Squared ($\mu_{13}$) & ', num2str(mu13,prec1), ' & (', num2str(mu13_SE,prec2), ') \\',...
    'Conditional - Father''s Education ($\mu_{14}$) & ', num2str(mu14,prec1), ' & (', num2str(mu14_SE,prec2), ') \\',...
    'Standard Deviation of Shock ($\sigma_{I}$) & ', num2str(sigma_NLI,prec1), ' & (', num2str(sigma_NLI_SE,prec2), ') \\ ']);

disp('\hline \hline');
disp('\end{tabular*}');
disp('\end{center}');
disp(['\textbf{Notes}: SEs are standard errors computed using a cluster bootstrap ',...
    'sampling each household with replacement.']);
disp('\end{table}');

diary off;

%% Graphs: plot the 3x3 Markov transition matrix probabilities

tvec = 1:fp.M;
markov_trans = obj_eval.markov_trans;

f=figure;
set(f,'visible','off');
for ii=1:3
    for jj=1:3
        subplot(3,3,jj+(ii-1)*3)
        plot(tvec, squeeze(markov_trans(ii,jj,1,:)),'k',tvec,squeeze(markov_trans(ii,jj,2,:)),'k--')
        xlabel('Age');
        ylabel('Prob.');
        axis([min(tvec) max(tvec) 0 1]);
    end
end
saveas(f,'figure_markov_all', 'png')
saveas(f,'figure_markov_all', 'eps')

%% Graph: Visualize initial conditions distributions for beta_{c,t0} at t0 = child_age97 when we start simulating the model

betac_init_draws = obj_eval.betac_init;

% Calculate average initial beta by age and parental education
tmpp = 1:(fp.M+1):size(data,1);
age1 = [3 6 8 10];
age2 = [5 7 9 11];
betac_init_graph = nan(length(age1),3); % rows = 3 age cat, columns = 3 educ cat: HS or less / SC or coll / graduat
for ii = 1:length(age1)
    tmpL = data(tmpp,col_long.child_age97)>=age1(ii) & data(tmpp,col_long.child_age97)<=age2(ii) & data(tmpp,col_long.dad_educ) <= 12;
    tmpM = data(tmpp,col_long.child_age97)>=age1(ii) & data(tmpp,col_long.child_age97)<=age2(ii) & data(tmpp,col_long.dad_educ) >= 13 & data(tmpp,col_long.dad_educ) <= 16;
    tmpH = data(tmpp,col_long.child_age97)>=age1(ii) & data(tmpp,col_long.child_age97)<=age2(ii) & data(tmpp,col_long.dad_educ) >= 17;

    tmp1 = betac_init_draws(tmpL,:,2);
    tmp2 = betac_init_draws(tmpM,:,2);
    tmp3 = betac_init_draws(tmpH,:,2);

    betac_init_graph(ii,1) = nanmean(tmp1(:));
    betac_init_graph(ii,2) = nanmean(tmp2(:));
    betac_init_graph(ii,3) = nanmean(tmp3(:));
end

f=figure;
set(f,'visible','off');
plot(age1,betac_init_graph(:,1),'k-',age1,betac_init_graph(:,2),'k-.',age1,betac_init_graph(:,3),'k--');
xlabel('Child'' Age in 1997');
ylabel('Average initial discount factor');
axis([min(age1) max(age1) 0 1]);
set(gca,'XTick',age1)
h = legend('Parent educ. - High school or less','Parent educ. - Some College or College','Parent educ. - Graduate','Location','southoutside');
set(gca,'XTickLabel',char('3-5','6-7','8-9','10-11'))
set(h, 'interpreter','latex');
saveas(f,'figure_betac_init_byageeduc', 'png')
saveas(f,'figure_betac_init_byageeduc', 'eps')

%% Table: Average and St.Dev. of Parents and Child Discount Factors, by age and parental education, Data vs. Simulation

% Create table of moments in data vs. simulation

% Load Osaka DATA moments (mean/std/median) by education category (bottom row = all hh):
osaka_moments_data = fp.osaka_moments_data; % see data_load_osaka.m
% Load Osaka SIM moments (mean/std/median) by education category (bottom row = all hh):
betap_moments_sim = fp.betap_moments_sim; % see data_load_osaka.m

% Load Steinberg DATA moments (mean/std/median) by education category (bottom row = all hh):
% first 4 rows = all ages 10-17, next 4 = ages 10-12, next 4 = ages 13-17.
steinberg_moments_data = fp.steinberg_moments_data; % see data_load_steinberg.m

% Define age cutoffs for child discount factors
% Note: in the Steinberg data, we only observe ages 10-17, whereas in the model, we simulate for ages 3-17
agevec1 = [3 3 6 10 13];
agevec2 = [17 5 9 12 17];

% define 3 EDUC categories: high school or less / some college or college / graduate
% Note: for PSID/CDS data that means
educvec1 = [0 13 17]; % must have same dimensions as fp.N_educ_betap and fp.N_educ_betac!
educvec2 = [12 16 20];
% For Steinberg and Osaka data sets, respectively, these cutoffs would be different due to differing SES variable definitions!

% Calculate "SIMULATED" moments of beta_c draws by age and education category (ALL, hs, sccoll, grad):
betac_moments_sim = NaN(length(age1)*(fp.N_educ_betac+1), 2); % rows = by age and SES cat / all hh, columns = mean/stdev

for ii = 1:length(agevec1)
    for jj = 1:length(educvec1)
        tmprows = sim_data(:,col_long.dad_educ) >= educvec1(jj) & sim_data(:,col_long.dad_educ) <= educvec2(jj) & sim_data(:,col_long.age) >= agevec1(ii) & sim_data(:,col_long.age) <= agevec2(ii);
        tmpc = sim_data(tmprows,col_long.n); % note: this doesn't load betas for age 17! only up to 16
        tmpbetas = tmpc./(1+tmpc);
        betac_moments_sim((ii-1)*(length(educvec1)+1)+jj,:) = [nanmean(tmpbetas) nanstd(tmpbetas)];
    end
    % now for all hh within the age category (averaging across all educ)
    tmprows = sim_data(:,col_long.age) >= agevec1(ii) & sim_data(:,col_long.age) <= agevec2(ii);
    tmpc = sim_data(tmprows,col_long.n); % note: this doesn't load betas for age 17! only up to 16
    tmpbetas = tmpc./(1+tmpc);
    betac_moments_sim((ii-1)*(length(educvec1)+1)+length(educvec1)+1,:) = [nanmean(tmpbetas) nanstd(tmpbetas)];
end
% % Verify output:
% osaka_moments_data
% betap_moments_sim
% 
% steinberg_moments_data
% betac_moments_sim

cd(paper_dir)
prec1 = '%.3f'; % precision for estimates

diary off;
diary_dir = [cd,'\table_fit_betas2.txt'];
force_delete(diary_dir);
diary('table_fit_betas2.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Discount Factors of Parents and Children - Data vs. Model Simulation}']);
disp(['\label{table:fit_betas}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcccc}']);
disp(['\hline \hline']);

disp(['\multicolumn{1}{l}{\textbf{(a) Parents/Adults}}  \rule{0pt}{3ex}& \multicolumn{2}{c}{\textbf{Data (Osaka-PPS)}} & \multicolumn{2}{c}{\textbf{Model Sim.}}\\']);
disp([' & \multicolumn{1}{c}{\textbf{Mean}} & \multicolumn{1}{c}{\textbf{St.Dev.}} & \multicolumn{1}{c}{\textbf{Mean}} & \multicolumn{1}{c}{\textbf{St.Dev.}}\\']);

disp(['All individuals \rule{0pt}{3ex} & ',...
    num2str(osaka_moments_data(end,1),prec1), ' & ', num2str(osaka_moments_data(end,2),prec1), ' & ',...
    num2str(betap_moments_sim(end,1),prec1), ' & ', num2str(betap_moments_sim(end,2),prec1), ' \\ ',...
    'High school or less & ',...
    num2str(osaka_moments_data(1,1),prec1), ' & ', num2str(osaka_moments_data(1,2),prec1), ' & ',...
    num2str(betap_moments_sim(1,1),prec1), ' & ', num2str(betap_moments_sim(1,2),prec1), ' \\ ',...
    'Some college or College & ',...
    num2str(osaka_moments_data(2,1),prec1), ' & ', num2str(osaka_moments_data(2,2),prec1), ' & ',...
    num2str(betap_moments_sim(2,1),prec1), ' & ', num2str(betap_moments_sim(2,2),prec1), ' \\ ',...
    'Graduate & ',...
    num2str(osaka_moments_data(3,1),prec1), ' & ', num2str(osaka_moments_data(3,2),prec1), ' & ',...
    num2str(betap_moments_sim(3,1),prec1), ' & ', num2str(betap_moments_sim(3,2),prec1), ' \\ ']);

disp(['\multicolumn{1}{l}{\textbf{(b) Children}}  \rule{0pt}{4ex}& \multicolumn{2}{c}{\textbf{Data (Steinberg et al.)}} & \multicolumn{2}{c}{\textbf{Model Sim.}}\\']);
disp([' & \multicolumn{1}{c}{\textbf{Mean}} & \multicolumn{1}{c}{\textbf{St.Dev.}} & \multicolumn{1}{c}{\textbf{Mean}} & \multicolumn{1}{c}{\textbf{St.Dev.}}\\']);

disp(['Ages 3-5 \rule{0pt}{3ex} & - & - & ', num2str(betac_moments_sim(8,1),prec1), ' & ', num2str(betac_moments_sim(8,2),prec1), ' \\ ',...
    '\quad High school or less & - & - & ',        num2str(betac_moments_sim(5,1),prec1), ' & ', num2str(betac_moments_sim(5,2),prec1), ' \\ ',...
    '\quad Some college or College & - & - & ',    num2str(betac_moments_sim(6,1),prec1), ' & ', num2str(betac_moments_sim(6,2),prec1), ' \\ ',...
    '\quad Graduate & - & - & ',                   num2str(betac_moments_sim(7,1),prec1), ' & ', num2str(betac_moments_sim(7,2),prec1), ' \\ ']);
disp(['Ages 6-9 \rule{0pt}{3ex} & - & - & ', num2str(betac_moments_sim(12,1),prec1), ' & ', num2str(betac_moments_sim(12,2),prec1), ' \\ ',...
    '\quad High school or less & - & - & ',        num2str(betac_moments_sim(9,1),prec1), ' & ', num2str(betac_moments_sim(9,2),prec1), ' \\ ',...
    '\quad Some college or College & - & - & ',    num2str(betac_moments_sim(10,1),prec1), ' & ', num2str(betac_moments_sim(10,2),prec1), ' \\ ',...
    '\quad Graduate & - & - & ',                   num2str(betac_moments_sim(11,1),prec1), ' & ', num2str(betac_moments_sim(11,2),prec1), ' \\ ']);
disp(['Ages 10-12 \rule{0pt}{3ex} & ',...
    num2str(steinberg_moments_data(8,1),prec1), ' & ', num2str(steinberg_moments_data(8,2),prec1), ' & ',...
    num2str(betac_moments_sim(16,1),prec1), ' & ', num2str(betac_moments_sim(16,2),prec1), ' \\ ',...
    '\quad High school or less & ',...
    num2str(steinberg_moments_data(5,1),prec1), ' & ', num2str(steinberg_moments_data(5,2),prec1), ' & ',...
    num2str(betac_moments_sim(13,1),prec1), ' & ', num2str(betac_moments_sim(13,2),prec1), ' \\ ',...
    '\quad Some college or College & ',...
    num2str(steinberg_moments_data(6,1),prec1), ' & ', num2str(steinberg_moments_data(6,2),prec1), ' & ',...
    num2str(betac_moments_sim(14,1),prec1), ' & ', num2str(betac_moments_sim(14,2),prec1), ' \\ ',...
    '\quad Graduate & ',...
    num2str(steinberg_moments_data(7,1),prec1), ' & ', num2str(steinberg_moments_data(7,2),prec1), ' & ',...
    num2str(betac_moments_sim(15,1),prec1), ' & ', num2str(betac_moments_sim(15,2),prec1), ' \\ ']);
disp(['Ages 13-17 \rule{0pt}{3ex} & ',...
    num2str(steinberg_moments_data(12,1),prec1), ' & ', num2str(steinberg_moments_data(12,2),prec1), ' & ',...
    num2str(betac_moments_sim(20,1),prec1), ' & ', num2str(betac_moments_sim(20,2),prec1), ' \\ ',...
    '\quad High school or less & ',...
    num2str(steinberg_moments_data(9,1),prec1), ' & ', num2str(steinberg_moments_data(9,2),prec1), ' & ',...
    num2str(betac_moments_sim(17,1),prec1), ' & ', num2str(betac_moments_sim(17,2),prec1), ' \\ ',...
    '\quad Some college or College & ',...
    num2str(steinberg_moments_data(10,1),prec1), ' & ', num2str(steinberg_moments_data(10,2),prec1), ' & ',...
    num2str(betac_moments_sim(18,1),prec1), ' & ', num2str(betac_moments_sim(18,2),prec1), ' \\ ',...
    '\quad Graduate & ',...
    num2str(steinberg_moments_data(11,1),prec1), ' & ', num2str(steinberg_moments_data(11,2),prec1), ' & ',...
    num2str(betac_moments_sim(19,1),prec1), ' & ', num2str(betac_moments_sim(19,2),prec1), ' \\ ']);  

disp('\hline \hline');
disp('\end{tabular*}');
disp('\end{center}');
disp(['\textbf{Notes}: Panel (a) shows the mean and standard deviation of adults'' annual discount factors, conditional on their own educational attainment. ',...
    'Panel (b) shows the mean and standard deviation of children''s annual discount factors, conditional on their age and their parents'' educational attainment. ',...
    'The data for adults stems from the 2010 U.S. survey of the Preference Parameters Study of Osaka University, for the subsample of adults between ages 25-65 (N = 4625). ',...
    'The data for children stems from the experimental surveys collected by Steinberg et al. (2009), where we restrict attention to children between ages 10-17 (N = 344). ',...
    'The first two columns show data moments, whereas the last two columns show summary statistics of the simulated (imputed) annual discount factors ',...
    'for our sample of PSID-CDS parents and children, respectively.']);

disp('\end{table}');

diary off;

%% Table: Sample fit of wages and non-labor income

prec = '%.2f'; % precision for moments

% Report key wage/NLI moments in the paper (Appendix has the full list)

% data moments for mothers wages
tmp_data_w1(1) = nanmean(data(:,col_long.mom_wage)); % uncond mean
tmp_data_w1(2) = nanstd(data(:,col_long.mom_wage)); % uncond std

% simulated moments for mothers wages
tmp_sim_w1(1) = nanmean(sim_data(:,col_long.mom_wage)); % uncond mean
tmp_sim_w1(2) = nanstd(sim_data(:,col_long.mom_wage)); % uncond std

% data moments for fathers wages
tmp_data_w2(1) = nanmean(data(:,col_long.dad_wage)); % uncond mean
tmp_data_w2(2) = nanstd(data(:,col_long.dad_wage)); % uncond std

% simulated moments for fathers wages
tmp_sim_w2(1) = nanmean(sim_data(:,col_long.dad_wage)); % uncond mean
tmp_sim_w2(2) = nanstd(sim_data(:,col_long.dad_wage)); % uncond std

% finally, select some relevant moments: see e.g. NLI_moment_funct.m

NLI_data = data(:,col_long.nonlabor_income); % observed data
NLI_sim = fp.NLI_sim; % simulated nonlabor income data

% discard some data
%data
dad_age_data = data(:,col_long.dad_age);
dad_age_data(dad_age_data<25 | dad_age_data>60) = NaN; % discard all men younger than 25 or older than 60 (too few observations)
dad_educ_data = data(:,col_long.dad_educ);
dad_educ_data(dad_educ_data<12) = NaN;
index_discard_data = (dad_age_data>0 & dad_educ_data>0);
NLI_data(index_discard_data == 0) = NaN; % discard these NLI observations
%sim
dad_age_sim = sim_data(:,col_long.dad_age);
dad_age_sim(dad_age_sim<25 | dad_age_sim>60) = NaN; % discard all men younger than 25 or older than 60 (too few observations)
dad_educ_sim = sim_data(:,col_long.dad_educ);
dad_educ_sim(dad_educ_sim<12) = NaN;
index_discard_sim = (dad_age_sim>0 & dad_educ_sim>0);
NLI_sim(index_discard_sim == 0) = NaN; % discard these NLI observations

% now differentiate between zeros and positives
% data
dummy_pos_data = NaN(size(NLI_data,1),1);
tmp0_data = find(NLI_data==0);
tmp1_data = find(NLI_data>0 & not(isnan(NLI_data)));
dummy_pos_data(tmp0_data) = 0;
dummy_pos_data(tmp1_data) = 1;
dummy_pos_data(index_discard_data == 0) = NaN;
% Note: this dummy is a NaN if NLI = NaN.

% sim
dummy_pos_sim = NaN(size(NLI_sim,1),1);
tmp0_sim = find(NLI_sim==0);
tmp1_sim = find(NLI_sim>0 & not(isnan(NLI_sim)));
dummy_pos_sim(tmp0_sim) = 0;
dummy_pos_sim(tmp1_sim) = 1;
dummy_pos_sim(index_discard_sim == 0) = NaN;
% Note: this dummy is a NaN if NLI = NaN.

% NLI data moments
tmp_data_NLI(1) = nanmean(NLI_data); % uncond mean
tmp_data_NLI(2) = nanstd(NLI_data); % uncond std
tmp_data_NLI(3) = nanmean(dummy_pos_data); % proportion of zeros

% NLI sim moments
tmp_sim_NLI(1) = nanmean(NLI_sim); % uncond mean
tmp_sim_NLI(2) = nanstd(NLI_sim); % uncond std
tmp_sim_NLI(3) = nanmean(dummy_pos_sim); % proportion of zeros

diary off;
diary_dir = [cd,'\table_fit_wages.txt'];
force_delete(diary_dir);
diary('table_fit_wages.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Sample Fit of Wages and Non-Labor Income}']);
disp(['\label{table:fit_wages}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcccc}']);
disp(['\hline \hline']);
disp(['\multicolumn{1}{l}{\textbf{(a) Hourly Wages}} \rule{0pt}{3ex} & \multicolumn{2}{c}{\textbf{Mother}} & \multicolumn{2}{c}{\textbf{Father}} \\']);
disp(['\cmidrule(l){2-3}\cmidrule(l){4-5}']);
disp(['& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}} & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}\\']);
disp(['\cline{2-5}']);

disp(['Average \rule{0pt}{3ex} & ', num2str(tmp_data_w1(1),prec), ' & ', num2str(tmp_sim_w1(1),prec), ' & ', num2str(tmp_data_w2(1),prec), ' & ', num2str(tmp_sim_w2(1),prec),' \\ ',...
    'Standard deviation &', num2str(tmp_data_w1(2),prec), ' & ', num2str(tmp_sim_w1(2),prec), ' & ', num2str(tmp_data_w2(2),prec), ' & ', num2str(tmp_sim_w2(2),prec), ' \\ ']);

disp(['\multicolumn{5}{l}{\textbf{(b) Weekly Non-labor Income}} \rule{0pt}{4ex}\\']);
disp(['& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}} & & \\']);
disp(['\cline{2-3}']);

disp(['Average \rule{0pt}{3ex} & ', num2str(tmp_data_NLI(1),prec), ' & ', num2str(tmp_sim_NLI(1),prec), ' \\',...
    'Standard deviation & ', num2str(tmp_data_NLI(2),prec), ' & ', num2str(tmp_sim_NLI(2),prec), ' \\',...
    'Fraction with $I>0$ & ', num2str(tmp_data_NLI(3),prec), ' & ', num2str(tmp_sim_NLI(3),prec), ' \\ ']);

disp('\hline \hline');
disp('\end{tabular*}');
disp('\end{center}');
disp(['\textbf{Notes}: Data is actual data. Simulated is the model prediction at estimated parameters.\\',...
    '\textbf{Source}: PSID-CDS combined sample from 1997, 2002 and 2007 interviews and PSID core data between 1986 and 2010.']);
disp('\end{table}');
diary off;

%% Table: Sample fit for time allocations

prec1 = '%.3f'; % precision for  LFP
prec = '%.2f'; % precision for rest

% define child age categories
beg_age = [3 6 9 13];
end_age = [5 8 12 16];

% a) prob employed
tmp_a_data = NaN*ones(4,2); % 4 age groups, 2 parents
tmp_a_sim = NaN*ones(4,2); % 4 age groups, 2 parents
% b) work hours if positive
tmp_b_data = NaN*ones(4,2); % 4 age groups, 2 parents
tmp_b_sim = NaN*ones(4,2); % 4 age groups, 2 parents
% c) active time with child
tmp_c_data = NaN*ones(4,2); % 4 age groups, 2 parents
tmp_c_sim = NaN*ones(4,2); % 4 age groups, 2 parents
% d) child self-investment time
tmp_d_data = NaN*ones(4,1); % 4 age groups
tmp_d_sim = NaN*ones(4,1); % 4 age groups

for t = 1:4
    % select the rows we need for this child age group
    row_ind_data = (data(:,col_long.age) >= beg_age(t) & data(:,col_long.age) <= end_age(t));
    row_ind_sim = (sim_data(:,col_long.age) >= beg_age(t) & sim_data(:,col_long.age) <= end_age(t));
    
    %a) prob employed
    % a.1. mothers
    temp1_data = nansum(data( row_ind_data,col_long.mom_hours) > 0);
    temp2_data = nansum((data( row_ind_data,col_long.mom_hours)) == 0);
    tmp_a_data(t,1) = temp1_data./(temp1_data + temp2_data);
    
    temp1_sim = nansum(sim_data( row_ind_sim,col_long.mom_hours) > 0);
    temp2_sim = nansum((sim_data( row_ind_sim,col_long.mom_hours)) == 0);
    tmp_a_sim(t,1) = temp1_sim./(temp1_sim + temp2_sim);
    
    % a.2. fathers
    temp1_data = nansum(data( row_ind_data,col_long.dad_hours) > 0);
    temp2_data = nansum((data( row_ind_data,col_long.dad_hours)) == 0);
    tmp_a_data(t,2) = temp1_data./(temp1_data + temp2_data);
    
    temp1_sim = nansum(sim_data( row_ind_sim,col_long.dad_hours) > 0);
    temp2_sim = nansum((sim_data( row_ind_sim,col_long.dad_hours)) == 0);
    tmp_a_sim(t,2) = temp1_sim./(temp1_sim + temp2_sim);
    
    %b) work hours if positive
    % b.1. mothers
    temp_data = (row_ind_data == 1 & data(:,col_long.mom_hours) > 0); % conditional on supplying positive hours
    temp_sim = (row_ind_sim == 1 & sim_data(:,col_long.mom_hours) > 0); % conditional on supplying positive hours
    tmp_b_data(t,1) = nanmean(data(temp_data,col_long.mom_hours));
    tmp_b_sim(t,1) = nanmean(sim_data(temp_sim,col_long.mom_hours));
    
    % b.2. fathers
    temp_data = (row_ind_data == 1 & data(:,col_long.dad_hours) > 0); % conditional on supplying positive hours
    temp_sim = (row_ind_sim == 1 & sim_data(:,col_long.dad_hours) > 0); % conditional on supplying positive hours
    tmp_b_data(t,2) = nanmean(data(temp_data,col_long.dad_hours));
    tmp_b_sim(t,2) = nanmean(sim_data(temp_sim,col_long.dad_hours));
    
    %c) active time with child
    % c.1. mothers
    tmp_c_data(t,1) = nanmean(data(row_ind_data,col_long.tau1));
    tmp_c_sim(t,1) = nanmean(sim_data(row_ind_sim,col_long.tau1));
    
    % c.2. fathers
    tmp_c_data(t,2) = nanmean(data(row_ind_data,col_long.tau2));
    tmp_c_sim(t,2) = nanmean(sim_data(row_ind_sim,col_long.tau2));
    
    % d) child active time
    tmp_d_data(t,1) = nanmean(data(row_ind_data,col_long.tau_c));
    tmp_d_sim(t,1) = nanmean(sim_data(row_ind_sim,col_long.tau_c));
end

diary off;
diary_dir = [cd,'\table_fit_time.txt'];
force_delete(diary_dir);
diary('table_fit_time.txt');
diary on;

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Sample Fit of Time Allocations by Child Age}']);
disp(['\label{table:fit_time}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} lcccc}']);
disp(['\hline \hline']);

disp('\multicolumn{5}{l}{\textbf{(a) Probability Work $>$ 0 Hours}} \rule{0pt}{3ex}\\');
disp([' & \multicolumn{2}{c}{\textbf{Mother}} & \multicolumn{2}{c}{\textbf{Father}}\\']);
disp('\cmidrule(l){2-3}\cmidrule(l){4-5}');
disp(['\textbf{Child Age} & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}'...
    '& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}\\']);
disp('\cline{1-5}');
disp(['3-5 \rule{0pt}{3ex}& ', num2str(tmp_a_data(1,1),prec1), ' & ', num2str(tmp_a_sim(1,1),prec1),' & ', num2str(tmp_a_data(1,2),prec1), ' & ', num2str(tmp_a_sim(1,2),prec1),' \\',...
    '6-8 & ', num2str(tmp_a_data(2,1),prec1), ' & ', num2str(tmp_a_sim(2,1),prec1),' & ', num2str(tmp_a_data(2,2),prec1), ' & ', num2str(tmp_a_sim(2,2),prec1),' \\',...
    '9-12 & ', num2str(tmp_a_data(3,1),prec1), ' & ', num2str(tmp_a_sim(3,1),prec1),' & ', num2str(tmp_a_data(3,2),prec1), ' & ', num2str(tmp_a_sim(3,2),prec1),' \\',...
    '13-16 & ', num2str(tmp_a_data(4,1),prec1), ' & ', num2str(tmp_a_sim(4,1),prec1),' & ', num2str(tmp_a_data(4,2),prec1), ' & ', num2str(tmp_a_sim(4,2),prec1),' \\ ']);

disp(['\hline \hline']);
disp('\multicolumn{5}{l}{\textbf{(b) Hours Worked if Work $>$ 0 Hours (Avg.)}} \rule{0pt}{4ex}\\');
disp([' & \multicolumn{2}{c}{\textbf{Mother}} & \multicolumn{2}{c}{\textbf{Father}}\\']);
disp('\cmidrule(l){2-3}\cmidrule(l){4-5}');
disp(['\textbf{Child Age} & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}'...
    '& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}\\']);
disp('\cline{1-5}');
disp(['3-5 \rule{0pt}{3ex}& ', num2str(tmp_b_data(1,1),prec), ' & ', num2str(tmp_b_sim(1,1),prec),' & ', num2str(tmp_b_data(1,2),prec), ' & ', num2str(tmp_b_sim(1,2),prec),' \\',...
    '6-8 & ', num2str(tmp_b_data(2,1),prec), ' & ', num2str(tmp_b_sim(2,1),prec),' & ', num2str(tmp_b_data(2,2),prec), ' & ', num2str(tmp_b_sim(2,2),prec),' \\',...
    '9-12 & ', num2str(tmp_b_data(3,1),prec), ' & ', num2str(tmp_b_sim(3,1),prec),' & ', num2str(tmp_b_data(3,2),prec), ' & ', num2str(tmp_b_sim(3,2),prec),' \\',...
    '13-16 & ', num2str(tmp_b_data(4,1),prec), ' & ', num2str(tmp_b_sim(4,1),prec),' & ', num2str(tmp_b_data(4,2),prec), ' & ', num2str(tmp_b_sim(4,2),prec),' \\ ']);

disp(['\hline \hline']);
disp('\multicolumn{5}{l}{\textbf{(c) Active Time with Child (Avg.)}} \rule{0pt}{4ex}\\');
disp([' & \multicolumn{2}{c}{\textbf{Mother}} & \multicolumn{2}{c}{\textbf{Father}}\\']);
disp('\cmidrule(l){2-3}\cmidrule(l){4-5}');
disp(['\textbf{Child Age} & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}'...
    '& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}}\\']);
disp('\cline{1-5}');
disp(['3-5 \rule{0pt}{3ex}& ', num2str(tmp_c_data(1,1),prec), ' & ', num2str(tmp_c_sim(1,1),prec),' & ', num2str(tmp_c_data(1,2),prec), ' & ', num2str(tmp_c_sim(1,2),prec),' \\',...
    '6-8 & ', num2str(tmp_c_data(2,1),prec), ' & ', num2str(tmp_c_sim(2,1),prec),' & ', num2str(tmp_c_data(2,2),prec), ' & ', num2str(tmp_c_sim(2,2),prec),' \\',...
    '9-12 & ', num2str(tmp_c_data(3,1),prec), ' & ', num2str(tmp_c_sim(3,1),prec),' & ', num2str(tmp_c_data(3,2),prec), ' & ', num2str(tmp_c_sim(3,2),prec),' \\',...
    '13-16 & ', num2str(tmp_c_data(4,1),prec), ' & ', num2str(tmp_c_sim(4,1),prec),' & ', num2str(tmp_c_data(4,2),prec), ' & ', num2str(tmp_c_sim(4,2),prec),' \\ ']);

disp(['\hline \hline']);
disp('\multicolumn{5}{l}{\textbf{(d) Child Self-Investment Time (Avg.)}} \rule{0pt}{4ex}\\');
disp(['\textbf{Child Age} \rule{0pt}{3ex}& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Simulated}} & & \\']);
disp('\cline{1-3}');
disp(['3-5 \rule{0pt}{3ex}& ', num2str(tmp_d_data(1,1),prec), ' & ', num2str(tmp_d_sim(1,1),prec),' & & \\',...
    '6-8 & ', num2str(tmp_d_data(2,1),prec), ' & ', num2str(tmp_d_sim(2,1),prec),' & & \\',...
    '9-12 & ', num2str(tmp_d_data(3,1),prec), ' & ', num2str(tmp_d_sim(3,1),prec),' & & \\',...
    '13-16 & ', num2str(tmp_d_data(4,1),prec), ' & ', num2str(tmp_d_sim(4,1),prec),' & & \\ ']);

disp('\hline \hline');
disp('\end{tabular*}');
disp('\end{center}');
disp(['\textbf{Notes}: Data is actual data. Simulated is the model prediction at estimated parameters.\\']);
disp(['\textbf{Source}: PSID-CDS combined sample from 1997, 2002 and 2007 interviews and PSID core data between 1986 and 2010.']);
disp('\end{table}');
diary off;

%% Figure: Average child's letter word raw score: data and model fit

% First prepare the data points we want to plot
xvec = [0:16]'; % kid age
yvec = NaN*ones(length(xvec),2);
for i = 1:size(xvec,1)
    age_i = xvec(i);
    % observed data
    rows = data(:,col_long.age)==age_i;
    yvec(i,1) = nanmean(data(rows,col_long.lw_raw));
    % sim data
    rows2 = sim_data(:,col_long.age)==age_i;
    yvec(i,2) = nanmean(sim_data(rows2,col_long.lw_raw));
end

% Data
f = figure;
set(f,'visible','off');
plot(xvec, yvec(:,1),'Color','k')
xlabel('Child''s age');
ylabel('Letter Word Score');
saveas(f,'figure_LWraw_data', 'png')
saveas(f,'figure_LWraw_data', 'eps')

% Sample fit
f = figure;
set(f,'visible','off');
plot(xvec, yvec(:,1),'-kd',xvec, yvec(:,2),'--k');
xlabel('Child''s age');
ylabel('Letter Word Score');
h = legend('Data','Simulation','Location','southoutside');
saveas(f,'figure_LWraw_fit', 'png')
saveas(f,'figure_LWraw_fit', 'eps')

%% Figures: Cognitive Technology parameters by child's age

xvec = [1:fp.M]'; % child age levels
yvec = NaN*ones(length(xvec),6);

% mom time
delta1_intercept = param_vector(40);
delta1_slope = param_vector(41);
delta1_educ = param_vector(84);
% dad time
delta2_intercept = param_vector(42);
delta2_slope = param_vector(43);
delta2_educ = param_vector(85);
% child time
delta3_intercept = param_vector(44);
delta3_slope = param_vector(45);
% child expenditures
delta4_intercept = param_vector(46);
delta4_slope = param_vector(47);
% lagged skills
delta5_intercept = param_vector(48);
delta5_slope = param_vector(49);
% TFP at each age
TFP_vec = [1; 1; param_vector(90:103)]; % ages1-2: TFP = 1, age 3-16: 14 parameters

delta1mat = NaN(fp.M,fp.N);
delta2mat = NaN(fp.M,fp.N);
for q = 1:fp.N
    data_here = data(1+(fp.M+1)*(q-1):(fp.M+1)*q, : );
    mom_educ = data_here(1,col_long.mom_educ);
    dad_educ = data_here(1,col_long.dad_educ);
    for t = 1:fp.M
        delta1mat(t,q) = exp(delta1_intercept + delta1_slope*t + delta1_educ*mom_educ);
        delta2mat(t,q) = exp(delta2_intercept + delta2_slope*t + delta2_educ*dad_educ);
    end % t
end % q

for t = 1:length(xvec)
    yvec(t,1) = nanmean(delta1mat(t,:));
    yvec(t,2) = nanmean(delta2mat(t,:));
    yvec(t,3) = exp(delta3_intercept + t * delta3_slope);
    yvec(t,4) = exp(delta4_intercept + t * delta4_slope);
    yvec(t,5) = exp(delta5_intercept + t * delta5_slope);
end

% Parental time productivity by education level
zvec = NaN*ones(length(xvec),8);
tlow = 12;
thigh = 16;
for t = 1:length(xvec)
    % for high school graduates
    zvec(t,1) = exp(delta1_intercept + t*delta1_slope + delta1_educ*tlow); % mother productivity
    zvec(t,2) = exp(delta2_intercept + t*delta2_slope + delta2_educ*tlow); % father productivity
    % for college graduates
    zvec(t,3) = exp(delta1_intercept + t*delta1_slope + delta1_educ*thigh); % mother productivity
    zvec(t,4) = exp(delta2_intercept + t*delta2_slope + delta2_educ*thigh); % father productivity
end

% Plot Technology parameters
% Time inputs
f = figure;
set(f,'visible','off');
plot(xvec(3:end), yvec(3:end,1),'-k',xvec(3:end), yvec(3:end,2),'--k',xvec(3:end), yvec(3:end,3),'-ko');
xlabel('Child''s Age');
ylabel('Input Productivity of Time');
h = legend('Mother''s Active Time $(\tau_1)$','Father''s Active Time $(\tau_2)$','Child''s Self-Investment Time $(\tau_c)$','Location','northeast');
set(h, 'interpreter','latex');
saveas(f,'figure_tech_time', 'png')
saveas(f,'figure_tech_time', 'eps')

% Parental time input productivity, for high school graduates vs college graduates
f = figure;
set(f,'visible','off');
plot(xvec(3:end), zvec(3:end,1),'-k+',xvec(3:end), zvec(3:end,3),'-kx',xvec(3:end), zvec(3:end,2),'-ko',xvec(3:end), zvec(3:end,4),'-kd');
h = legend('Mother''s Active Time $(\tau_1)$ - High School','Mother''s Active Time $(\tau_1)$ - College','Father''s Active Time $(\tau_2)$ - High School','Father''s Active Time $(\tau_2)$ - College','Location','northeast');
set(h, 'interpreter','latex');
ylabel('Input Productivity of Time');
saveas(f,'figure_tech_time_byeduc', 'png')
saveas(f,'figure_tech_time_byeduc', 'eps')

% Child goods input
f = figure;
set(f,'visible','off');
plot(xvec(3:end), yvec(3:end,4),'-k');
xlabel('Child''s Age');
ylabel('Input Productivity');
set(h, 'interpreter','latex');
saveas(f,'figure_tech_childgoods', 'png')
saveas(f,'figure_tech_childgoods', 'eps')

% Previous child quality
f = figure;
set(f,'visible','off');
plot(xvec(3:end), yvec(3:end,5),'-k');
xlabel('Child''s Age');
ylabel('Input Productivity');
set(h, 'interpreter','latex');
saveas(f,'figure_tech_laggedq', 'png')
saveas(f,'figure_tech_laggedq', 'eps')

% Alternatively: plot log(TFP) instead
f = figure;
set(f,'visible','off');
plot(xvec(3:end), log(TFP_vec(3:end)),'-k');
xlabel('Child''s Age');
ylabel('Log(TFP)');
saveas(f,'figure_logTFP', 'png')
saveas(f,'figure_logTFP', 'eps')

%% Make some extra graphs to show the model fit (for Appendix)

graphs_data_sim(fp,col_long,data,sim_data,paper_dir)

%% Table 8: detailed summary of reported school times

school_tot = data(:,col_long.tau_school); % total reported school times in CDS (cat8+9; regular and other)
% school_reg = data(:,col_long.tau_school_reg); % regular school time
% school_oth = data(:,col_long.tau_school_oth); % other school time
% tmp = [school_tot school_reg school_oth];
tmp = school_tot;

captions = {'Total School Time $s_t$ by Child Age'};

tmp2 = 3:fp.M;
table = NaN(length(tmp2),9);

school=tmp;
for j = 1:length(tmp2)
    k=tmp2(j);
    rows_j = data(:,col_long.age) == k;
    var_j = school(rows_j);
    table(j,1) = nanmean(var_j);
    table(j,2) = nanstd(var_j);
    table(j,3) = nanmin(var_j);
    table(j,4) = prctile(var_j,25);
    table(j,5) = prctile(var_j,50);
    table(j,6) = prctile(var_j,75);
    table(j,7) = nanmax(var_j);
    table(j,8) = sum(var_j == 0); % nr of zeros
    table(j,9) = sum(not(isnan(var_j))); % nr of observations
end
% disp('Table with summary statistics of school time (copy-pasted in paper Appendix):')
% table

%% Create formatted tables of data vs simulated moments

tmp = [data_mom(fp.index_moments_1) sim_mom(fp.index_moments_1)]; % 68 moments, all used (time inputs, labor supply, wages)
tmp2 = [data_mom(fp.index_moments_2) sim_mom(fp.index_moments_2)]; % 14 moments, all used (LW by each age)
tmp3 = [data_mom(fp.index_moments_3) sim_mom(fp.index_moments_3)]; % 15 moments, all used, wages by age and educ, 2nd order wage moments
tmp4 = [data_mom(fp.index_moments_LW) sim_mom(fp.index_moments_LW)]; % 14 moments, all used, mean and stdev of LW and autocorrelations
tmp5 = [data_mom(fp.index_moments_tau1) sim_mom(fp.index_moments_tau1)]; % 4 moments, corr of tau1 and LW
tmp6 = [data_mom(fp.index_moments_tau2) sim_mom(fp.index_moments_tau2)]; % 4 moments, corr of tau2 and LW
tmp7 = [data_mom(fp.index_moments_tauc) sim_mom(fp.index_moments_tauc)]; % 4 moments, corr of tauc and LW
tmp8 = [data_mom(fp.index_moments_educ) sim_mom(fp.index_moments_educ)]; % 12 moments, corr of LW/educ1/educ2/tau1/tau2/tauc
tmp9 = [data_mom(fp.index_moments_CCT) sim_mom(fp.index_moments_CCT)]; % 7 moments, corr of CCT and other things
tmp10 = [data_mom(fp.index_moments_external) sim_mom(fp.index_moments_external)]; % 10 moments, all used, external expenditure/beta_c moments
% total: 152 moments used in estimation

tmp11 = [fp.NLI_data_mom fp.NLI_sim_mom]; % another 57 moments for external non-labor income process - estimated separately
% see NLI_moment_funct.m to see how the moments are ordered and defined

diary off;
diary_dir = [cd,'\table_moments_new1.txt'];
force_delete(diary_dir);
diary('table_moments_new1.txt');
diary on;

prec = '%.2f'; % precision for this table

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Data vs. Simulated Moments - Part 1}']);
disp(['\label{table:moments_new1}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lcclcc}']);
disp(['\hline \hline']);
disp([' & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Sim.}} & ',...
    '& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Sim.}} \\']);
disp(['\hline']);

disp(['$E(\tau_{1,t}|t \in [3,5])$ \rule{0pt}{3ex} & ',...
    num2str(tmp(17,1),prec),' & ',num2str(tmp(17,2),prec),' & ', ...
    '$\sigma(\tau_{1,t}|t \in [3,5])$ & ',...
    num2str(tmp(21,1),prec),' & ',num2str(tmp(21,2),prec),' \\ ', ...
    '$E(\tau_{1,t}|t \in [6,8])$ & ',...
    num2str(tmp(18,1),prec),' & ',num2str(tmp(18,2),prec),' & ', ...
    '$\sigma(\tau_{1,t}|t \in [6,8])$ & ',...
    num2str(tmp(22,1),prec),' & ',num2str(tmp(22,2),prec),' \\ ', ...
    '$E(\tau_{1,t}|t \in [9,12])$ &',...
    num2str(tmp(19,1),prec),' & ',num2str(tmp(19,2),prec),' & ', ...
    '$\sigma(\tau_{1,t}|t \in [9,12])$ &',...
    num2str(tmp(23,1),prec),' & ',num2str(tmp(23,2),prec),' \\ ', ...
    '$E(\tau_{1,t}|t \in [13,16])$ & ',...
    num2str(tmp(20,1),prec),' & ',num2str(tmp(20,2),prec),' & ', ...
    '$\sigma(\tau_{1,t}|t \in [13,16])$ & ',...
    num2str(tmp(24,1),prec),' & ',num2str(tmp(24,2),prec),' \\ ']);

disp(['$E(\tau_{2,t}|t \in [3,5])$ \rule{0pt}{3ex} & ',...
    num2str(tmp(25,1),prec),' & ',num2str(tmp(25,2),prec),' & ', ...
    '$\sigma(\tau_{2,t}|t \in [3,5])$ & ',...
    num2str(tmp(29,1),prec),' & ',num2str(tmp(29,2),prec),' \\ ', ...
    '$E(\tau_{2,t}|t \in [6,8])$ & ',...
    num2str(tmp(26,1),prec),' & ',num2str(tmp(26,2),prec),' & ', ...
    '$\sigma(\tau_{2,t}|t \in [6,8])$ & ',...
    num2str(tmp(30,1),prec),' & ',num2str(tmp(30,2),prec),' \\ ', ...
    '$E(\tau_{2,t}|t \in [9,12])$ &',...
    num2str(tmp(27,1),prec),' & ',num2str(tmp(27,2),prec),' & ', ...
    '$\sigma(\tau_{2,t}|t \in [9,12])$ &',...
    num2str(tmp(31,1),prec),' & ',num2str(tmp(31,2),prec),' \\ ', ...
    '$E(\tau_{2,t}|t \in [13,16])$ & ',...
    num2str(tmp(28,1),prec),' & ',num2str(tmp(28,2),prec),' & ', ...
    '$\sigma(\tau_{2,t}|t \in [13,16])$ & ',...
    num2str(tmp(32,1),prec),' & ',num2str(tmp(32,2),prec),' \\ ']);

disp(['$E(\tau_{c,t}|t \in [3,5])$ \rule{0pt}{3ex} & ',...
    num2str(tmp(33,1),prec),' & ',num2str(tmp(33,2),prec),' & ', ...
    '$\sigma(\tau_{c,t}|t \in [3,5])$ & ',...
    num2str(tmp(37,1),prec),' & ',num2str(tmp(37,2),prec),' \\ ', ...
    '$E(\tau_{c,t}|t \in [6,8])$ & ',...
    num2str(tmp(34,1),prec),' & ',num2str(tmp(34,2),prec),' & ', ...
    '$\sigma(\tau_{c,t}|t \in [6,8])$ & ',...
    num2str(tmp(38,1),prec),' & ',num2str(tmp(38,2),prec),' \\ ', ...
    '$E(\tau_{c,t}|t \in [9,12])$ &',...
    num2str(tmp(35,1),prec),' & ',num2str(tmp(35,2),prec),' & ', ...
    '$\sigma(\tau_{c,t}|t \in [9,12])$ &',...
    num2str(tmp(39,1),prec),' & ',num2str(tmp(39,2),prec),' \\ ', ...
    '$E(\tau_{c,t}|t \in [13,16])$ & ',...
    num2str(tmp(36,1),prec),' & ',num2str(tmp(36,2),prec),' & ', ...
    '$\sigma(\tau_{c,t}|t \in [13,16])$ & ',...
    num2str(tmp(40,1),prec),' & ',num2str(tmp(40,2),prec),' \\ ']);

disp(['$E(h_{1,t}|t \in [3,5],h_{1,t}>0)$ \rule{0pt}{3ex} & ',...
    num2str(tmp(1,1),prec),' & ',num2str(tmp(1,2),prec),' & ', ...
    '$\sigma(h_{1,t}|t \in [3,5],h_{1,t}>0)$ & ',...
    num2str(tmp(5,1),prec),' & ',num2str(tmp(5,2),prec),' \\ ', ...
    '$E(h_{1,t}|t \in [6,8],h_{1,t}>0)$ & ',...
    num2str(tmp(2,1),prec),' & ',num2str(tmp(2,2),prec),' & ', ...
    '$\sigma(h_{1,t}|t \in [6,8],h_{1,t}>0)$ & ',...
    num2str(tmp(6,1),prec),' & ',num2str(tmp(6,2),prec),' \\ ', ...
    '$E(h_{1,t}|t \in [9,12],h_{1,t}>0)$ &',...
    num2str(tmp(3,1),prec),' & ',num2str(tmp(3,2),prec),' & ', ...
    '$\sigma(h_{1,t}|t \in [9,12],h_{1,t}>0)$ &',...
    num2str(tmp(7,1),prec),' & ',num2str(tmp(7,2),prec),' \\ ', ...
    '$E(h_{1,t}|t \in [13,16],h_{1,t}>0)$ & ',...
    num2str(tmp(4,1),prec),' & ',num2str(tmp(4,2),prec),' & ', ...
    '$\sigma(h_{1,t}|t \in [13,16],h_{1,t}>0)$ & ',...
    num2str(tmp(8,1),prec),' & ',num2str(tmp(8,2),prec),' \\ ']);

disp(['$E(h_{2,t}|t \in [3,5],h_{2,t}>0)$ \rule{0pt}{3ex} & ',...
    num2str(tmp(9,1),prec),' & ',num2str(tmp(9,2),prec),' & ', ...
    '$\sigma(h_{2,t}|t \in [3,5],h_{2,t}>0)$ & ',...
    num2str(tmp(13,1),prec),' & ',num2str(tmp(13,2),prec),' \\ ', ...
    '$E(h_{2,t}|t \in [6,8],h_{2,t}>0)$ & ',...
    num2str(tmp(10,1),prec),' & ',num2str(tmp(10,2),prec),' & ', ...
    '$\sigma(h_{2,t}|t \in [6,8],h_{2,t}>0)$ & ',...
    num2str(tmp(14,1),prec),' & ',num2str(tmp(14,2),prec),' \\ ', ...
    '$E(h_{2,t}|t \in [9,12],h_{2,t}>0)$ & ',...
    num2str(tmp(11,1),prec),' & ',num2str(tmp(11,2),prec),' & ', ...
    '$\sigma(h_{2,t}|t \in [9,12],h_{2,t}>0)$ & ',...
    num2str(tmp(15,1),prec),' & ',num2str(tmp(15,2),prec),' \\ ', ...
    '$E(h_{2,t}|t \in [13,16],h_{2,t}>0)$ & ',...
    num2str(tmp(12,1),prec),' & ',num2str(tmp(12,2),prec),' & ', ...
    '$\sigma(h_{2,t}|t \in [13,16],h_{2,t}>0)$ & ',...
    num2str(tmp(16,1),prec),' & ',num2str(tmp(16,2),prec),' \\ ']);

disp(['$Pr(h_{1,t} > 0|t \in [3,5])$ \rule{0pt}{3ex} & ',...
    num2str(tmp(57,1),prec),' & ',num2str(tmp(57,2),prec),' & ', ...
    '$Pr(h_{2,t} > 0|t \in [3,5])$ & ',...
    num2str(tmp(61,1),prec),' & ',num2str(tmp(61,2),prec),' \\ ', ...
    '$Pr(h_{1,t} > 0|t \in [6,8])$ & ',...
    num2str(tmp(58,1),prec),' & ',num2str(tmp(58,2),prec),' & ', ...
    '$Pr(h_{2,t} > 0|t \in [6,8])$ & ',...
    num2str(tmp(62,1),prec),' & ',num2str(tmp(62,2),prec),' \\ ', ...
    '$Pr(h_{1,t} > 0|t \in [9,12])$ & ',...
    num2str(tmp(59,1),prec),' & ',num2str(tmp(59,2),prec),' & ', ...
    '$Pr(h_{2,t} > 0|t \in [9,12])$ & ',...
    num2str(tmp(63,1),prec),' & ',num2str(tmp(63,2),prec),' \\ ', ...
    '$Pr(h_{1,t} > 0|t \in [13,16])$ & ',...
    num2str(tmp(60,1),prec),' & ',num2str(tmp(60,2),prec),' & ', ...
    '$Pr(h_{2,t} > 0|t \in [13,16])$ & ',...
    num2str(tmp(64,1),prec),' & ',num2str(tmp(64,2),prec),' \\ ',...
    '$Pr((h_{1,t},h_{2,t})>0|t \in [3,5])$ & ',...
    num2str(tmp(65,1),prec),' & ',num2str(tmp(65,2),prec),' & ', ...
    '$Pr((h_{1,t},h_{2,t})>0|t \in [6,8])$ & ',...
    num2str(tmp(66,1),prec),' & ',num2str(tmp(66,2),prec),' \\ ', ...
    '$Pr((h_{1,t},h_{2,t})>0|t \in [9,12])$ & ',...
    num2str(tmp(67,1),prec),' & ',num2str(tmp(67,2),prec),' & ', ...
    '$Pr((h_{1,t},h_{2,t})>0|t \in [13,16])$ & ',...
    num2str(tmp(68,1),prec),' & ',num2str(tmp(68,2),prec),' \\ ']);

disp(['$E(w_{1,t}|age_{1,t} \in [18,32])$ \rule{0pt}{3ex} & ',...
    num2str(tmp(41,1),prec),' & ',num2str(tmp(41,2),prec),' & ', ...
    '$\sigma(w_{1,t}|age_{1,t} \in [18,32])$ & ',...
    num2str(tmp(45,1),prec),' & ',num2str(tmp(45,2),prec),' \\ ', ...
    '$E(w_{1,t}|age_{1,t} \in [33,37])$  & ',...
    num2str(tmp(42,1),prec),' & ',num2str(tmp(42,2),prec),' & ', ...
    '$\sigma(w_{1,t}|age_{1,t} \in [33,37])$ & ',...
    num2str(tmp(46,1),prec),' & ',num2str(tmp(46,2),prec),' \\ ', ...
    '$E(w_{1,t}|age_{1,t} \in [38,42])$  & ',...
    num2str(tmp(43,1),prec),' & ',num2str(tmp(43,2),prec),' & ', ...
    '$\sigma(w_{1,t}|age_{1,t} \in [38,42])$ & ',...
    num2str(tmp(47,1),prec),' & ',num2str(tmp(47,2),prec),' \\ ', ...
    '$E(w_{1,t}|age_{1,t} \in [43,65])$ & ',...
    num2str(tmp(44,1),prec),' & ',num2str(tmp(44,2),prec),' & ', ...
    '$\sigma(w_{1,t}|age_{1,t} \in [43,65])$ & ',...
    num2str(tmp(48,1),prec),' & ',num2str(tmp(48,2),prec),' \\ ']);

disp(['$E(w_{2,t}|age_{1,t} \in [18,32])$ \rule{0pt}{3ex} & ',...
    num2str(tmp(49,1),prec),' & ',num2str(tmp(49,2),prec),' & ', ...
    '$\sigma(w_{2,t}|age_{1,t} \in [18,32])$ & ',...
    num2str(tmp(53,1),prec),' & ',num2str(tmp(53,2),prec),' \\ ', ...
    '$E(w_{2,t}|age_{1,t} \in [33,37])$  & ',...
    num2str(tmp(50,1),prec),' & ',num2str(tmp(50,2),prec),' & ', ...
    '$\sigma(w_{2,t}|age_{1,t} \in [33,37])$ & ',...
    num2str(tmp(54,1),prec),' & ',num2str(tmp(54,2),prec),' \\ ', ...
    '$E(w_{2,t}|age_{1,t} \in [38,42])$  & ',...
    num2str(tmp(51,1),prec),' & ',num2str(tmp(51,2),prec),' & ', ...
    '$\sigma(w_{2,t}|age_{1,t} \in [38,42])$ & ',...
    num2str(tmp(55,1),prec),' & ',num2str(tmp(55,2),prec),' \\ ', ...
    '$E(w_{2,t}|age_{1,t} \in [43,65])$ & ',...
    num2str(tmp(52,1),prec),' & ',num2str(tmp(52,2),prec),' & ', ...
    '$\sigma(w_{2,t}|age_{1,t} \in [43,65])$ & ',...
    num2str(tmp(56,1),prec),' & ',num2str(tmp(56,2),prec),' \\ ']);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes:} The table shows moments related to productive time investments and labor supply conditional on the child''s age, ' ...
    'and hourly accepted wages conditional on the parent''s age.' ]);
disp(['\end{table}']);

diary off;

diary_dir = [cd,'\table_moments_new2.txt'];
force_delete(diary_dir);
diary('table_moments_new2.txt');
diary on;

prec = '%.2f'; % precision for this table

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Data vs. Simulated Moments - Part 2}']);
disp(['\label{table:moments_new2}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lcclcc}']);
disp(['\hline \hline']);
disp([' & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Sim.}} & ',...
    '& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Sim.}} \\']);
disp(['\hline']);

% other wage moments
disp(['$E(w_{1,t}|s_{1} = 12, age_{1,t} < 40)$ \rule{0pt}{3ex} & ',...
    num2str(tmp3(1,1),prec),' & ',num2str(tmp3(1,2),prec),' & ', ...
    '$E(w_{1,t}|s_{1} = 12, age_{1,t} \geq 40)$ & ',...
    num2str(tmp3(2,1),prec),' & ',num2str(tmp3(2,2),prec),' \\ ', ...
    '$E(w_{1,t}|s_{1} \in[13,15], age_{1,t} < 40)$ & ',...
    num2str(tmp3(3,1),prec),' & ',num2str(tmp3(3,2),prec),' & ', ...
    '$E(w_{1,t}|s_{1} \in[13,15], age_{1,t} \geq 40)$ & ',...
    num2str(tmp3(4,1),prec),' & ',num2str(tmp3(4,2),prec),' \\ ', ...
    '$E(w_{1,t}|s_{1} \geq 16, age_{1,t} < 40)$ & ',...
    num2str(tmp3(5,1),prec),' & ',num2str(tmp3(5,2),prec),' & ', ...
    '$E(w_{1,t}|s_{1} \geq 16, age_{1,t} \geq 40)$ & ',...
    num2str(tmp3(6,1),prec),' & ',num2str(tmp3(6,2),prec),' \\ ', ...
    '$E(w_{2,t}|s_{2} = 12, age_{2,t} < 40)$ & ',...
    num2str(tmp3(7,1),prec),' & ',num2str(tmp3(7,2),prec),' & ', ...
    '$E(w_{2,t}|s_{2} = 12, age_{2,t} \geq 40)$ & ',...
    num2str(tmp3(8,1),prec),' & ',num2str(tmp3(8,2),prec),' \\ ', ...
    '$E(w_{2,t}|s_{2} \in[13,15], age_{2,t} < 40)$ & ',...
    num2str(tmp3(9,1),prec),' & ',num2str(tmp3(9,2),prec),' & ', ...
    '$E(w_{2,t}|s_{2} \in[13,15], age_{2,t} \geq 40)$ & ',...
    num2str(tmp3(10,1),prec),' & ',num2str(tmp3(10,2),prec),' \\ ', ...
    '$E(w_{2,t}|s_{2} \geq 16, age_{2,t} < 40)$ & ',...
    num2str(tmp3(11,1),prec),' & ',num2str(tmp3(11,2),prec),' & ', ...
    '$E(w_{2,t}|s_{2} \geq 16, age_{2,t} \geq 40)$ & ',...
    num2str(tmp3(12,1),prec),' & ',num2str(tmp3(12,2),prec),' \\ ',...
    '$\sigma(w_{1,t})$ & ',...
    num2str(tmp3(13,1),prec),' & ',num2str(tmp3(13,2),prec),' & ', ...
    '$corr(w_{1,t},w_{2,t})$  & ',...
    num2str(tmp3(15,1),prec),' & ',num2str(tmp3(15,2),prec),' \\ ', ...
    '$\sigma(w_{2,t})$ & ',...
    num2str(tmp3(14,1),prec),' & ',num2str(tmp3(14,2),prec),' & & & \\ ']);

disp(['$E(LW_t | t = 3)$ \rule{0pt}{3ex} & ',...
    num2str(tmp2(1,1),prec),' & ',num2str(tmp2(1,2),prec),' & ', ...
    '$E(LW_t)$, All obs. & ',...
    num2str(tmp4(1,1),prec),' & ',num2str(tmp4(1,2),prec),' \\ ', ...
    '$E(LW_t | t = 4)$ & ',...
    num2str(tmp2(2,1),prec),' & ',num2str(tmp2(2,2),prec),' & ', ...
    '$E(LW_t| t \in [3,7])$ & ',...
    num2str(tmp4(2,1),prec),' & ',num2str(tmp4(2,2),prec),' \\ ', ...
    '$E(LW_t | t = 5)$ & ',...
    num2str(tmp2(3,1),prec),' & ',num2str(tmp2(3,2),prec),' & ', ...
    '$E(LW_t | t \in [8,11])$ & ',...
    num2str(tmp4(3,1),prec),' & ',num2str(tmp4(3,2),prec),' \\ ', ...
    '$E(LW_t | t = 6)$ & ',...
    num2str(tmp2(4,1),prec),' & ',num2str(tmp2(4,2),prec),' & ', ...
    '$E(LW_t| t \in [12,16])$ & ',...
    num2str(tmp4(4,1),prec),' & ',num2str(tmp4(4,2),prec),' \\ ', ...
    '$E(LW_t | t = 7)$ & ',...
    num2str(tmp2(5,1),prec),' & ',num2str(tmp2(5,2),prec),' & ', ...
    '$\sigma(LW_t)$, All obs. & ',...
    num2str(tmp4(5,1),prec),' & ',num2str(tmp4(5,2),prec),' \\ ', ...
    '$E(LW_t | t = 8)$ & ',...
    num2str(tmp2(6,1),prec),' & ',num2str(tmp2(6,2),prec),' & ', ...
    '$\sigma(LW_t| t \in [3,7])$ & ',...
    num2str(tmp4(6,1),prec),' & ',num2str(tmp4(6,2),prec),' \\ ', ...
    '$E(LW_t | t = 9)$ & ',...
    num2str(tmp2(7,1),prec),' & ',num2str(tmp2(7,2),prec),' & ', ...
    '$\sigma(LW_t| t \in [8,11])$ & ',...
    num2str(tmp4(7,1),prec),' & ',num2str(tmp4(7,2),prec),' \\ ', ...
    '$E(LW_t | t = 10)$ & ',...
    num2str(tmp2(8,1),prec),' & ',num2str(tmp2(8,2),prec),' & ', ...
    '$\sigma(LW_t| t \in [12,16])$ & ',...
    num2str(tmp4(8,1),prec),' & ',num2str(tmp4(8,2),prec),' \\ ', ...
    '$E(LW_t | t = 11)$ & ',...
    num2str(tmp2(9,1),prec),' & ',num2str(tmp2(9,2),prec),' & ', ...
    '$E(\Delta LW_{t,t+5})$, All obs. & ',...
    num2str(tmp4(9,1),prec),' & ',num2str(tmp4(9,2),prec),' \\ ', ...
    '$E(LW_t | t = 12)$ & ',...
    num2str(tmp2(10,1),prec),' & ',num2str(tmp2(10,2),prec),' & ', ...
    '$E(\Delta LW_{t,t+5} | t \in [3,7])$ & ',...
    num2str(tmp4(10,1),prec),' & ',num2str(tmp4(10,2),prec),' \\ ', ...
    '$E(LW_t | t = 13)$ & ',...
    num2str(tmp2(11,1),prec),' & ',num2str(tmp2(11,2),prec),' & ', ...
    '$E(\Delta LW_{t,t+5} | t \in [8,11])$ & ',...
    num2str(tmp4(11,1),prec),' & ',num2str(tmp4(11,2),prec),' \\ ', ...
    '$E(LW_t | t = 14)$ & ',...
    num2str(tmp2(12,1),prec),' & ',num2str(tmp2(12,2),prec),' & ', ...
    '$corr(LW_t,LW_{t+5})$, All obs. & ',...
    num2str(tmp4(12,1),prec),' & ',num2str(tmp4(12,2),prec),' \\ ', ...
    '$E(LW_t | t = 15)$ & ',...
    num2str(tmp2(13,1),prec),' & ',num2str(tmp2(13,2),prec),' & ', ...
    '$corr(LW_t,LW_{t+5}| t \in [3,7])$ & ',...
    num2str(tmp4(13,1),prec),' & ',num2str(tmp4(13,2),prec),' \\ ', ...
    '$E(LW_t | t = 16)$ & ',...
    num2str(tmp2(14,1),prec),' & ',num2str(tmp2(14,2),prec),' & ', ...
    '$corr(LW_t,LW_{t+5}| t \in [8,11])$ & ',...
    num2str(tmp4(14,1),prec),' & ',num2str(tmp4(14,2),prec),' \\ ']);

disp(['$corr(\tau_{1,t},LW_t)$, All ages \rule{0pt}{3ex} & ',...
    num2str(tmp5(1,1),prec),' & ',num2str(tmp5(1,2),prec),' & ', ...
    '$corr(\tau_{1,t},\Delta LW_{t,t+5})$, All ages & ',...
    num2str(tmp5(2,1),prec),' & ',num2str(tmp5(2,2),prec),' \\ ', ...
    '$corr(\tau_{1,t},\Delta LW_t) | t \in [3,7])$ & ',...
    num2str(tmp5(3,1),prec),' & ',num2str(tmp5(3,2),prec),' & ', ...
    '$corr(\tau_{1,t},\Delta LW_t) | t \in [8,11])$ & ',...    
    num2str(tmp5(4,1),prec),' & ',num2str(tmp5(4,2),prec),' \\ ', ...
    '$corr(\tau_{2,t},LW_t)$, All ages \rule{0pt}{3ex} & ',...
    num2str(tmp6(1,1),prec),' & ',num2str(tmp6(1,2),prec),' & ', ...
    '$corr(\tau_{2,t},\Delta LW_{t,t+5})$, All ages & ',...
    num2str(tmp6(2,1),prec),' & ',num2str(tmp6(2,2),prec),' \\ ', ...
    '$corr(\tau_{2,t},\Delta LW_t) | t \in [3,7])$ & ',...
    num2str(tmp6(3,1),prec),' & ',num2str(tmp6(3,2),prec),' & ', ...
    '$corr(\tau_{2,t},\Delta LW_t) | t \in [8,11])$ & ',...    
    num2str(tmp6(4,1),prec),' & ',num2str(tmp6(4,2),prec),' \\ ', ...
    '$corr(\tau_{c,t},LW_t)$, All ages \rule{0pt}{3ex} & ',...
    num2str(tmp7(1,1),prec),' & ',num2str(tmp7(1,2),prec),' & ', ...
    '$corr(\tau_{c,t},\Delta LW_{t,t+5})$, All ages & ',...
    num2str(tmp7(2,1),prec),' & ',num2str(tmp7(2,2),prec),' \\ ', ...
    '$corr(\tau_{c,t},\Delta LW_t) | t \in [3,7])$ & ',...
    num2str(tmp7(3,1),prec),' & ',num2str(tmp7(3,2),prec),' & ', ...
    '$corr(\tau_{c,t},\Delta LW_t) | t \in [8,11])$ & ',...    
    num2str(tmp7(4,1),prec),' & ',num2str(tmp7(4,2),prec),' \\ ']);

% correlations between educ/LW and educ/taus
disp(['$corr(s_1,LW_t|t \in [3,8])$ \rule{0pt}{3ex} & ',...
    num2str(tmp8(1,1),prec),' & ',num2str(tmp8(1,2),prec),' & ', ...
    '$corr(s_2,LW_t|t \in [3,8])$  & ',...
    num2str(tmp8(3,1),prec),' & ',num2str(tmp8(3,2),prec),' \\ ', ...
    '$corr(s_1,LW_t|t \in [9,12])$  & ',...
    num2str(tmp8(5,1),prec),' & ',num2str(tmp8(5,2),prec),' & ', ...
    '$corr(s_2,LW_t|t \in [9,12])$  & ',...
    num2str(tmp8(7,1),prec),' & ',num2str(tmp8(7,2),prec),' \\ ', ...
    '$corr(s_1,LW_t|t \in [13,16])$  & ',...
    num2str(tmp8(9,1),prec),' & ',num2str(tmp8(9,2),prec),' & ', ...
    '$corr(s_2,LW_t|t \in [13,16])$  & ',...
    num2str(tmp8(11,1),prec),' & ',num2str(tmp8(11,2),prec),' \\ ', ...
    '$corr(s_1,\tau_{1,t}|t \in [3,8])$  \rule{0pt}{3ex} & ',...
    num2str(tmp8(2,1),prec),' & ',num2str(tmp8(2,2),prec),' & ', ...
    '$corr(s_2,\tau_{2,t}|t \in [3,8])$  & ',...
    num2str(tmp8(4,1),prec),' & ',num2str(tmp8(4,2),prec),' \\ ', ...
    '$corr(s_1,\tau_{1,t}|t \in [9,12])$  & ',...
    num2str(tmp8(6,1),prec),' & ',num2str(tmp8(6,2),prec),' & ', ...
    '$corr(s_2,\tau_{2,t}|t \in [3,8])$  & ',...
    num2str(tmp8(8,1),prec),' & ',num2str(tmp8(8,2),prec),' \\ ', ...
    '$corr(s_1,\tau_{1,t}|t \in [13,16])$  & ',...
    num2str(tmp8(10,1),prec),' & ',num2str(tmp8(10,2),prec),' & ', ...
    '$corr(s_2,\tau_{2,t}|t \in [13,16])$  & ',...
    num2str(tmp8(12,1),prec),' & ',num2str(tmp8(12,2),prec),' \\ ']);
    
disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes:} The table shows moments related to hourly accepted wages conditional on parental age and schooling (where $12$ years denotes high school, ' ...
    '$[13,15]$ is some college and $16+$ is college or more), Letter Word scores (i.e., levels, ',...
    'changes over time and autocorrelations conditional on the child''s age), ', ...
    'correlations between time investments and (changes in) test scores; ',...
    'between parental schooling and test scores; and between parental schooling and time investments.']);
disp(['\end{table}']);

diary off;

diary_dir = [cd,'\table_moments_new3.txt'];
force_delete(diary_dir);
diary('table_moments_new3.txt');
diary on;

prec = '%.2f'; % precision for this table

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Data vs. Simulated Moments - Part 3}']);
disp(['\label{table:moments_new3}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lcclcc}']);
disp(['\hline \hline']);
disp([' & \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Sim.}} & ',...
    '& \multicolumn{1}{c}{\textbf{Data}} & \multicolumn{1}{c}{\textbf{Sim.}} \\']);
disp(['\hline']);

% CCT related moments
disp(['$E(CCT|t \in [8,16])$ \rule{0pt}{3ex} & ',...
    num2str(tmp9(1,1),prec),' & ',num2str(tmp9(1,2),prec),' & ', ...
    '$corr(CCT_t,age_t|t \in [8,16])$  & ',...
    num2str(tmp9(2,1),prec),' & ',num2str(tmp9(2,2),prec),' \\ ', ...
    '$corr(CCT_t,s_1|t \in [8,16])$  & ',...
    num2str(tmp9(3,1),prec),' & ',num2str(tmp9(3,2),prec),' & ', ...
    '$corr(CCT_t,s_2|t \in [8,16])$  & ',...
    num2str(tmp9(4,1),prec),' & ',num2str(tmp9(4,2),prec),' \\ ', ...
    '$corr(CCT_t,\tau_{1,t}|t \in [8,16])$  & ',...
    num2str(tmp9(5,1),prec),' & ',num2str(tmp9(5,2),prec),' & ', ...
    '$corr(CCT_t,\tau_{2,t}|t \in [8,16])$ & ',...
    num2str(tmp9(6,1),prec),' & ',num2str(tmp9(6,2),prec),' \\ ', ...
    '$corr(CCT_t,LW_t|t \in [8,16])$  & ',...
    num2str(tmp9(7,1),prec),' & ',num2str(tmp9(7,2),prec),' & & & \\ ']);

% external moments on child expenditures and child betas
disp(['$E(e_t+x_t | Y_t \in (1000,2500])$ \rule{0pt}{3ex} & ',...
    num2str(tmp10(1,1),prec),' & ',num2str(tmp10(1,2),prec),' & ', ...
    '$E(\frac{e_t+x_t}{Y_t}|Y_t \in (1000,2500])$  & ',...
    num2str(tmp10(2,1),prec),' & ',num2str(tmp10(2,2),prec),' \\ ', ...
    '$E(\beta_{c,t}|t \in [10,12])$ \rule{0pt}{3ex} & ',...
    num2str(tmp10(3,1),prec),' & ',num2str(tmp10(3,2),prec),' & ', ...
    '$E(\beta_{c,t}|t \in [13,17])$ & ',...
    num2str(tmp10(7,1),prec),' & ',num2str(tmp10(7,2),prec),' \\ ', ...
    '$E(\beta_{c,t}|t \in [10,12], s_2 \leq 12)$  & ',...
    num2str(tmp10(4,1),prec),' & ',num2str(tmp10(4,2),prec),' & ', ...
    '$E(\beta_{c,t}|t \in [13,17], s_2 \leq 12)$ & ',...
    num2str(tmp10(8,1),prec),' & ',num2str(tmp10(8,2),prec),' \\ ', ...
    '$E(\beta_{c,t}|t \in [10,12], s_2 \in [13,16])$  & ',...
    num2str(tmp10(5,1),prec),' & ',num2str(tmp10(5,2),prec),' & ', ...
    '$E(\beta_{c,t}|t \in [13,17], s_2 \in [13,16])$ & ',...
    num2str(tmp10(9,1),prec),' & ',num2str(tmp10(9,2),prec),' \\ ', ...
    '$E(\beta_{c,t}|t \in [10,12], s_2 \geq 17)$  & ',...
    num2str(tmp10(6,1),prec),' & ',num2str(tmp10(6,2),prec),' & ', ...
    '$E(\beta_{c,t}|t \in [13,17], s_2 \geq 17)$ & ',...
    num2str(tmp10(10,1),prec),' & ',num2str(tmp10(10,2),prec),' \\ ']);

% moments of Nonlabor income process
disp(['$E(I_t)$ \rule{0pt}{3ex} & ',...
    num2str(tmp11(15,1),prec),' & ',num2str(tmp11(15,2),prec),' & ', ...
    '$\sigma(I_t)$  & ',...
    num2str(tmp11(16,1),prec),' & ',num2str(tmp11(16,2),prec),' \\ ', ...
    '$E(I_t|I_t>0)$ & ',...
    num2str(tmp11(1,1),prec),' & ',num2str(tmp11(1,2),prec),' & ', ...
    '$\sigma(I_t|I_t>0)$  & ',...
    num2str(tmp11(2,1),prec),' & ',num2str(tmp11(2,2),prec),' \\ ', ...
    '$E(\mathbbm{1}[I_t>0])$ & ',...
    num2str(tmp11(8,1),prec),' & ',num2str(tmp11(8,2),prec),' & ', ...
    '$\sigma(\mathbbm{1}[I_t>0])$  & ',...
    num2str(tmp11(9,1),prec),' & ',num2str(tmp11(9,2),prec),' \\ ', ...
    '$E(I_t | age_{1,t} \in [25,32])$  \rule{0pt}{3ex} & ',...
    num2str(tmp11(46,1),prec),' & ',num2str(tmp11(46,2),prec),' & ', ...
    '$E(I_t | age_{2,t} \in [25,33])$  & ',...
    num2str(tmp11(22,1),prec),' & ',num2str(tmp11(22,2),prec),' \\ ', ...
    '$E(I_t | age_{1,t} \in [33,37])$ & ',...
    num2str(tmp11(47,1),prec),' & ',num2str(tmp11(47,2),prec),' & ', ...
    '$E(I_t | age_{2,t} \in [34,38])$  & ',...
    num2str(tmp11(23,1),prec),' & ',num2str(tmp11(23,2),prec),' \\ ', ...
    '$E(I_t | age_{1,t} \in [38,42])$ & ',...
    num2str(tmp11(48,1),prec),' & ',num2str(tmp11(48,2),prec),' & ', ...
    '$E(I_t | age_{2,t} \in [39,43])$  & ',...
    num2str(tmp11(24,1),prec),' & ',num2str(tmp11(24,2),prec),' \\ ', ...
    '$E(I_t | age_{1,t} \in [43,55] )$ & ',...
    num2str(tmp11(49,1),prec),' & ',num2str(tmp11(49,2),prec),' & ', ...
    '$E(I_t | age_{2,t} \in [44,60])$  & ',...
    num2str(tmp11(25,1),prec),' & ',num2str(tmp11(25,2),prec),' \\ ', ...
    '$E(I_t | I_t>0, age_{1,t} \in [25,32])$  \rule{0pt}{3ex} & ',...
    num2str(tmp11(50,1),prec),' & ',num2str(tmp11(50,2),prec),' & ', ...
    '$E(I_t | I_t>0, age_{2,t} \in [25,33])$  & ',...
    num2str(tmp11(26,1),prec),' & ',num2str(tmp11(26,2),prec),' \\ ', ...
    '$E(I_t | I_t>0, age_{1,t} \in [33,37])$ & ',...
    num2str(tmp11(51,1),prec),' & ',num2str(tmp11(51,2),prec),' & ', ...
    '$E(I_t | I_t>0, age_{2,t} \in [34,38])$  & ',...
    num2str(tmp11(27,1),prec),' & ',num2str(tmp11(27,2),prec),' \\ ', ...
    '$E(I_t | I_t>0, age_{1,t} \in [38,42])$ & ',...
    num2str(tmp11(52,1),prec),' & ',num2str(tmp11(52,2),prec),' & ', ...
    '$E(I_t | I_t>0, age_{2,t} \in [39,43])$  & ',...
    num2str(tmp11(28,1),prec),' & ',num2str(tmp11(28,2),prec),' \\ ', ...
    '$E(I_t | I_t>0, age_{1,t} \in [43,55] )$ & ',...
    num2str(tmp11(53,1),prec),' & ',num2str(tmp11(53,2),prec),' & ', ...
    '$E(I_t | I_t>0, age_{2,t} \in [44,60])$  & ',...
    num2str(tmp11(29,1),prec),' & ',num2str(tmp11(29,2),prec),' \\ ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{1,t} \in [25,32])$  \rule{0pt}{3ex} & ',...
    num2str(tmp11(54,1),prec),' & ',num2str(tmp11(54,2),prec),' & ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{2,t} \in [25,33])$  & ',...
    num2str(tmp11(30,1),prec),' & ',num2str(tmp11(30,2),prec),' \\ ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{1,t} \in [33,37])$ & ',...
    num2str(tmp11(55,1),prec),' & ',num2str(tmp11(55,2),prec),' & ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{2,t} \in [34,38])$  & ',...
    num2str(tmp11(31,1),prec),' & ',num2str(tmp11(31,2),prec),' \\ ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{1,t} \in [38,42])$ & ',...
    num2str(tmp11(56,1),prec),' & ',num2str(tmp11(56,2),prec),' & ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{2,t} \in [39,43])$  & ',...
    num2str(tmp11(32,1),prec),' & ',num2str(tmp11(32,2),prec),' \\ ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{1,t} \in [43,55] )$ & ',...
    num2str(tmp11(57,1),prec),' & ',num2str(tmp11(57,2),prec),' & ', ...
    '$E(\mathbbm{1}[I_t>0] | age_{2,t} \in [44,60])$  & ',...
    num2str(tmp11(33,1),prec),' & ',num2str(tmp11(33,2),prec),' \\ ']);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes:} The table shows moments related to the use of CCTs and household characteristics or choices; expenditures on children; ',...
    'the distribution of children''s discount factors conditional on child age and parental education; ',...
    'and the exogenous non-labor income process (a subset of moments pertaining to the nonlabor income process, ',...
    'including various interactions terms between $I_t$ and parental characteristics, are omitted for brevity).']);
    
disp(['\end{table}']);

diary off;
